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Set Books 


G. D. Smith, Numerical Solution of Partial Differential Equations (Oxford. 1971). 
H. F. Weinberger, A First Course in Partial Differential Equations (Xerox. 1965). 


It is essential to have these books; the course is based on them and will not make 
sense without them. They are referred to in the text as S and W respectively. 


Unit 11 is based on S: Chapter 5, pages 131 to 137, 143 to 149. 


Conventions 


Before working through this text make sure you have read A Guide to the Course: 
Partial Differential Equations of Applied Mathematics. References to Open University 
courses in mathematics take the form: 


Unit M100 13, Integration I for the Mathematics Foundation Course. 


Unit M20f 23, The Wave Equation for the Linear Mathematics Course. 


11.0 INTRODUCTION 


In Unit 5, Initial Value Problems and Unit 8, Stability we dealt with the finite-difference 
method for solving initial value problems in one space dimension. We now turn 
our attention to pure boundary value problems in two space dimensions. Here we 
seek the solution of a partial differential equation in some finite region of the xy-plane 
bounded by a closed curve on which a (boundary) condition is specified at every 
point. À typical problem on which we shall concentrate throughout this unit is 
Poisson's equation, 

eu ou 

=a tz Fl: 

6x ey 


in some rectangular region. 


We begin by applying the finite-difference techniques of Unit 5 to two problems which 
illustrate that there are two distinct topics for investigation in the numerical solution 
of boundary value problems. The first and most obvious task is to search for good 
numerical methods for solving the finite-difference equations which arise. The second 
task is familiar from our work in Unit 8, namely the investigation of the accuracy of 
the finite-difference solution as a solution to the boundary value problem. This 
latter investigation will, as always, raise the questions of stability and convergence. 
In Section 11.2 we shall look at these questions in some detail for Poisson's equation, 
where we shall find that a maximum principle applies just as in the analysis of boundary 
value problems for differential equations in Unit 3, Elliptic and Parabolic Equations. 
The maximum principle will enable us to find some bounds on the error in our 
numerical solution. 


The remainder of the unit deals with the numerical solution of the finite-diflerence 
equations. Our two examples highlight the important fact that a finite-difference 
method applied to boundary value problems gives rise to very large sets of simul- 
taneous equations which are linear if the differential equation is linear. The associated 
matrices are of a special sparse form, which means that they contain a significant 
majority of zero elements. We know from Unit M201 8, Numerical Solution of 
Simultaneous Algebraic Equations, that linear equations can be solved by various 
direct methods, such as Gauss elimination. With such methods it is sometimes possible 
to take sparsity into account, as for example in Unit 5 where we modified the elimina- 
tion method to give a recurrence relation for solving tridiagonal systems of equations. 
Unfortunately the sparse matrices associated with elliptic partial differential equations 
have a more complicated structure, and it is difficult to find economic direct methods 
of solution. 


There is, however, another important class of methods for solving linear equations, 
namely iterative methods which we first considered in Unit M /00 28, Linear Algebra 
IV. It turns out that iterative methods are specially suitable for our particular class 
of large sparse matrices, and we shall look at iterative processes in general terms as 
well as discussing specific methods. 
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11.1 TWO EXAMPLES 


The first reading passage begins with a review of some general properties of boundary 


value problems which were discussed in detail in W: Sections 10 and 11 


Notes 
(i) 


(ii) 


(iii) 


(iv) 


(vi) 


S: page 132, line 9 

A specified é@/én on C does not determine ¢ uniquely in S. If $ is a solution, 
then $ + A is a solution for any constant A (see W: page 54). 

S: page 133, line 5 

Although we are now dealing with pure boundary value problems we still 
approximate the differential equation by finite-difference schemes similar to 
those for initial value problems derived in Unit 5. Employing the usual central- 
difference replacement for second-order derivatives we can write the required 
finite-difference scheme as 


ET. 4536 42-0, 


where we use the same mesh spacing h in both the x- and y-directions. 


S: page 134, line 1 
In order to improve the accuracy we could increase the number of mesh points 
(thereby reducing h). This increases the number of simultaneous equations which 
have to be solved. 
S: page 134, line 3 
It can be shown that the solution of Poisson's equation near a re-entrant corner 
(where the interior angle exceeds 180°) is oscillatory, making our finite-difference 
methods inaccurate there. (See L. Fox, ed.. Numerical Solution of Ordinary and 
Partial Differential Equations, p. 303.) 
S: page 134, lines 12 to 17 
You need not worry about the derivation of the equation on line 15, but note that 
putting y = 45y (provided 4 is a positive constant) yields 

Ou du q 

Ox Oy K 


which is Poisson's equation. 


S: page 134, lines —8 and —7 
The symmetry condition yields 
u=0 onx ==1 
and 
ðu ii 
Ta u ony = —1. 


Tt also implies that there will be no heat flow across Oy and Ox; hence, 


Qu 
mas 0 on Oy 
Ox 

and 
[27] 
4-0 on Ox. 
oy 


(vii) S: page 136, lines — 2 and —1 


The matrix A is an example of a type of matrix which occurs frequently in this 
subject. In its partitioned form we see that A has a tridiagonal appearance and 


we therefore refer to it as a block tridiagonal matrix. Note also that each of the 
diagonal blocks in A is a tridiagonal matrix. 
SAQ I 


If, in Equation (5.2) on S: page 134, we employ the transformation y- BF to 
obtain Poisson's equation (note (v) above), what changes must be made to the bound- 
ary values to retain the same problem? What differences will there be in the numerical 
solution of this new problem compared to the solution given in S: pages 134 to 137, 
if the same mesh is used? 


(Solution on p. 35.) 


SAQ 2 
S: page 161, Exercise 1 (Omit the last part of the question which is on page 162.) 


(Solution on p. 35.) 


SAQ3 


S: page 162, Exercise 2 
Note that the mesh point numbering given in the question is different from that in 
Fig. 5.2; the reference to Fig. 5.2 is for illustration only. 


(Solution on p. 36.) 
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11.2 THE FIVE-POINT FORMULA 


11.2.0 Introduction 


The two examples of the previous section have illustrated two important facts 
associated with the solution of boundary value problems by finite-difference methods. 
The first is that the boundary value problem is approximated by an algebraic problem 
in which it is necessary to solve a system of simultaneous linear equations. The 
second fact is that the system of equations is generally very large. Therefore, our 
analysis of the numerical process divides into two distinct parts. We first need to know 
the conditions under which the solution of the approximating algebraic problem 
approaches the true solution of the differential equation; and then we need to investi- 
gate efficient methods for solving the algebraic problem. S deals only with solving the 
algebraic problem. This section, therefore, deals with the question of how well the 
solution of the algebraic problem approximates the true solution of the boundary 
value problem. There is no single technique which answers this in general. We have 
therefore chosen to investigate in some detail the application of one simple finite- 
difference scheme to the solution of Poisson's equation, which can be dealt with using 
a maximum principle much like the one we studied in Unit 3. 


We look at the problem 


a2 2 

U E «f duel (la) 
ey? 

U=g (xyeC, (1b) 


where D is some domain in the xy-plane and C is its boundary. 


We cover D by a rectangular mesh such that the mesh spacing in the x-direction is 
Ax and in the y-direction is Ay. We shall approximate the Laplacian by central 
differences to give the five-point formula, 


SS os wy litiy 7 Dui + Minty | Majer — Dg tija, 
ax? Oy ny (Ax)? ap o 
Hence the finite-difference problem can be written as 
[Li]; = fi; ije Ds, (2a) 
ug — 8g LEC, (2b) 


. where D, is the set of mesh points interior to D, C, is the set of mesh points* on C 
and L is the finite-diflerence operator given by 


o2 ò? 
[Lu]; = las * ais LN 


IESE — 2uj; + Wi-j y Magen — Quay + Uy jay 


(Ax)? (Ay)? (3) 


Usually we take Ax = Ay = h (say) in which case the finite-difference formula reduces 
to 


t 
[Lu]; = qa Uag P tiaj t Uji + uis, — Aut; j). 


In Section 11.2.1 we show that the solution to the finite-difference problem (2) with 
L defined by Equation (3) is unique. We do this by establishing a maximum ii zi l 
for this problem which enables us to show that the homogeneous noob tert a s 
= gi; = 0, has only the trivial solution t; j = 0. In consequence the solutio he 
€ nonhomogeneous problem is unique. The maximum principle has a hes 


imporlant use in that we can determine from it a bound on our approximate soluti 
ximate s ion. 


We shall obtain this bound in Section 11.22, and we deduce that the chosen finit 
a e- 


difference scheme is convergent to the solution the differe 
olution of i ation 
Ax, Ay ~ 0. erential equation as 


* We shall assume throughout that C coincides with a union of mesh lines. 
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11.2.4 The Maximum Principle 


In Unit 3 (Section 3.2.1) we saw that if U is a solution of 


eu eu 
ex? oy? 


=f  (xyeD 


and f (x. y) > 0 in D, then if U is continuous on D ù C it attains its maximum on C. 
Further, in the case where f(x, y) = 0, which gives Laplace's equation, both the 
maximum and minimum values of U are attained on C. A desirable property of the 
solution of the finite-difference problem, if it is to be a reasonable approximation 
to the solution of the boundary value problem. is that it too satisfies a similar 
maximum principle. This is in fact the case as we show in the following theorem. 


THEOREM 1 
Let L be the finite-difference operator 
85 
(Ax — (Ay)? 
(a) If v; is a function defined on the set of mesh points Dy Y C, which satisfies 


Le 


4 BO for all i, je D4. 
then 
mpra, € mati 
(b) If v; satisfies 


Lv; SO for all ije Da, 


UM 
ij 
then 


min vj 2 _ tij. 


Proof 


(a) The proof is by contradiction. Assume that at some point Pa with coordinates 
(igAx, joAy) we have 


Viajo = M 
where 

M > tij forall ije D, 
and 

M > vij for all ije C4. 


That is to say, we assume that the maximum value M of the function v is attained 
at some point Py in the interior of the domain and that the boundary values are 
all less than M. 


We can write the five-point formula as 


Ua. l3 + Uy 3 I 1 
poe ee e deca 4 
Leo = a? t (AXE [ass * aye P 


where we have used the notation 


Co = ij 
Ci = Vint tint 
02 = Yig- liar 


Ua = Cig date 


Us = Vin ssa 


as in the following diagram. 


M-v9X 


(Ax)? + (Ay)? l 


from Equation (4). But M > v; j in D, implies, by Equation (5), that 
v. = M r= 1.2,3,4, 


since the coefficients of the v, are all positive and their sum is 1. (For if, say, 
v, < M then the right-hand side of (5) is less than M.) 


This means that if v9 € D, is a maximum then so are v; , v5, v, and v,. We now 
choose one of the v, (r = 1, 2, 3, 4) as a new central point for the five-point formula 
and repeat the above argument. For example, choosing v, asthe new central point 
will involve vo, vs, vg and v; in the five-point formula (see figure). Since we have 
just shown that v, is a maximum, so are vs, vg and vy. In this way we can 
include all points in D, and C, to show that if one point in the interior is a 
maximum then all points are maxima, that is 


vj; = M at all points of D, and C,. 
This contradicts the assumption that v « M on C, and part (a) is proved. 
(b) Anargument similar to (a) could be employed here but it is simpler to note that 
max(—v,) = —min v; 
and 
L[-v],; = — Luj. 


Hence, if v; ; satisfies the hypothesis of (b) that Lv, ; € 0 then — 9; satisfies the 
hypothesis of (a). The conclusion of part (a) for —v,; is identical to the con- 
clusion of (b) for bij. 


Hence we have established the maximum and minimum principles for the five-point 
finite-difference replacement of the Laplacian. 


If we now consider the homogeneous form of Equations (2) in which 7 
for all i je D, U C, then Theorem | tells us that both the maximum a 
values of v; ; are zero, and therefore the only solution is the trivial one. 


ij— 8i5=0 
nd minimum 


This fact enables us to prove uniqueness. For, Suppose that Equations ( 
solutions u and v, and that w = u — p. Then, 
Luj = fi, lje Ds, 


Mig = Big LjeCs, 


2) have two 


and 


il je Ds, 


i. je C,. 
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and therefore by subtraction 
Lw;,,; = 0 ijeD,, 
wyj=0 ijeCa. 


The above observation that the homogeneous case has only the trivial solution, 
w = 0, implies that u = v, and the solution of Equations (2) is therefore unique. 


11.2.2 Error Analysis 


When the finite-difference equations (2) are written down for every internal mesh 
point we have a set of simultaneous equations 


Au — b, (6) 


where the elements of the column vector u are the mesh values U; j, A is a square 
matrix of order equal to the number of internal mesh points and whose coefficients 
depend on the finite-difference operator L, and b is a column vector of known elements 
depending on the internal mesh values of f; ; and the boundary values g; ;. If U is the 
corresponding column vector of values of the true solution of the differential equation, 
then we define the local truncation error of our finite-difference approximation at 


the point i,j by the formula 
n = LU;,; = Sig. 


If we have no computer storage errors (rounding errors) in the numbers J; j, g; j and 
uj j it is then easy to see that the vector U satisfies, instead of (6), the equation 


AU=b+T ' (7) 


where T is the column vector of local truncation errors. The column vector e of 
global errors e; ; where 


eij = Uj tij 
is thus given by 
Ae — T, (8a) 


by subtracting Equation (6) from Equation (7). We have already shown that Equation 
(6) has a unique solution, so that the inverse matrix A~! exists, and we can write 
Equation (8a) in the form 


e= A`'T. (8b) 


Now, just as in initial value problems, we would hope that as we steadily reduced the 
mesh spacings Ax and Ay our computed solutions would get closer to the true 
solution. Both the matrix A and (as we shall show) the vector T depend on Ax and Ay, 
and what we are hoping to show is that lim A~ 'T = 0 (the zero vector) as Ax, Ay ~0*. 
This will certainly happen if lim T = 0 as Ax, Ay ~0*, and the inverse matrix A7 ' 
is bounded as Ax, Ay ~0* (in the sense discussed in Unit 8 Section 8.1). The first 
requirement, that lim T = 0 as Ax, Ay ~0*, means (as before) that our finite- 
difference approximations are consistent with the given differential equation. The 
requirement of a bounded A^! is somewhat analogous to the concept of stability 
in our treatment of initial value problems, and if lime = 0 as Ax. Ay ~0* we say 
that our method is convergent, again in analogy with the corresponding initial value 
situation. 


We can obtain the local truncation error T; ; of our five-point finite-difference formula 
applied to Equation (1a) by Taylor's Theorem as explained in Unit 5. We find (see 
SAQ 4) that 


i ] + O((Ax)*) + O((Ayy*). (9) 


It follows that if the various derivatives exist throughout the domain D then lim T;; «0 


as Ax and Ay ~0*, and we have consistency. 


It is also convenient to simplify expression (9) for T; ; and find an upper Bound for it. 
If we use Taylors Theorem with remainder (see Unit 8, note (ii) of Section 8.2.1) it is 
easy to show (see SAQ 4) that 


IT, jS BIM (Ax? + MAy]. 


where M, and M, are the largest values of |@*U/@x*| and J8*Uj8y^| in DUC. 


(10) 


Let us now turn our attention to the global error. As a consequence of the following 
theorem, which uses the result of Theorem 1, we shall obtain a bound on the global 
error ej j. 


THEOREM 2 


Let v be any function defined on the set of mesh points D, Y C, contained in some 
rectangular region (0, a] x (0. b]. Then 


max |v] € max |v| + 4a? max |Le], 
Ds Cy Da 
where L is the five-point finite-difference operator 


(Ax — (Ay)? 


VA 


Proof 

We introduce the function on D, o C, given by 
Qij = Max), 

and observe that for all i je D, U Ca, 


OS, S 4a? 


(11) 
and 
Loi = ilü + 1? — 2? + (i - 1) 
=1. 
Now define two functions v* and v^ on D, U C, by 
vt = +04 Nó, 
(12) 


where 


N= max |Lp; jl. 


Clearly, operating on both sides of Equation (12) by L, 
Lož, = tLe; +N >0 lje D, 


the inequality coming from the definition of N. Applying the maximum principle 
(Theorem 1a) to each of v* we obtain 


vx max vé 
a 


< max[+ v; No; by Equation (12) 
a 

< max[ v; ] + 4Na? by Equation (11). 
" 


The definition of vå and the fact that Pij 2 0 imply that 

to; vé, VijeD, U Ca- 
Hence, 

tvy S max[ toi] + iNa* 

^ 
& max lol + 3a? N. 

Since the right-hand side of the final inequality is independent of i,j in D,, the theorem 
follows. i 


[Note that we can replace $a? by $b? since the function y given by y; j = jAy}? 
can be used to replace ¢;,; in the proof. Thus we can reduce the bound by choosing 
the smaller of a and 5.] 
The bound on the error is now easily found. The local truncation error T; ; is given by 
LU,;— Jij 27, ije Ds, 
where U, ; is the true solution of the boundary value problem. Hence, 
Lej = LU, — Lu; = Ti; ije Ds, 
and also, 
eij = UU; —u-0 ijeC,, 


since the same boundary data are used in the calculation of the true solution U and 
the finite-difference solution u. Applying Theorem 2 to the function e = U — u, we 
obtain 


Lie 
max engl € ya max[T, j. (13) 


The consequence of Equation (13) is that, provided the 7;; are finite, the e;j are 
bounded. 


Finally, our bound (10) for the truncation error shows that 
maxle, | < 3a" [M«(Ax)* + M,(Ay)?) ' (14) 


and if Ax = Ay = h, and M is the greater of M, and M,, we have the elegant result 


maxle; | < 3;«^h* M, (15) 


showing that the global error is of the same order as the local truncation error. Since 
the right-hand side of (15) approaches zero as h ^» 0, our numerical scheme is con- 
vergent. 


It is, of course, difficult to compute this error bound, since it depends on the maximum 
value of the fourth derivatives of the function we are trying to compute! In general 
we do not try to do this, but use (15) to guarantee convergence and to indicate the 
rate of convergence. For example, if we halve the interval h we shall reduce the error 
to approximately one quarter of its previous value. and this is important and useful 
information. 


In all this analysis we have assumed that there are no computer errors in the stored 
values of the elements of A and b in the linear equation (6). that is, of the coeflicients 


a 5 s E 
of the operator L and the quantities f; , and gij- Possible uncertainties 5 Serg 
errors in f; ; and g;j can be analysed by using the maximum principle. If Ẹ is the 
computed solution of the problem then we can write 


Li = fij Rij i. je D, 
ty = Bij + Rij LjeCa. 
where we have distinguished between the errors R' made at points in D, and those 
R" made in approximating the boundary data. 
As before, we can obtain the expressions 
L(U-iü);- T; Rij  ijeD, 
U;;—-i;- -Riy injec, 
and applying Theorem 2 to the function U — ii, we obtain 
max|U -ü <s max|R"] + ta max|T — R!. 
and once again the error is bounded. 


Alternatively we could argue that possible small uncertainties in the data will not 
cause significant inherent instability since our problem is properly posed, as we 
have seen in Unit 3. 


We have also assumed that we introduce very little induced instability in solving 
the linear equations. We already know (Unit M201 8) that this is true if we use the 
direct method of Gauss elimination with pivoting. It can also be shown, though we 
shall not prove this, that no additional induced instability need be associated with 
our preferred methods of iteration (for our particular boundary value problems): we 
shall consider these methods in Section 11.3. 


SAQ 4 
Derive an expression for the local truncation error of the five-point formula applied to 


SU FU 
ax? ay? =, ' 


in a bounded domain D, and obtain a bound on its absolute value. 


(Solution on p. 37.) 


11.3 GENERAL ITERATIVE METHODS 


11.3.1 Discussion 


There are two classes of numerical methods for solving systems of linear algebraic 
equations, direct and iterative. In direct methods, of which Gaussian elimination 
with pivoting is the main example, the solution is obtained in a known number of 
arithmetic operations (which depends on the order of the matrix). With exact arith- 
metic the solution so obtained is exact, and only the unavoidable occurrence of digital 
computer rounding errors impairs this. In iterative methods, on the other hand, we 
obtain solutions by a process of successive approximation. We shall rarely obtain an 
exact solution, even with exact arithmetic, but if the method converges then we shall 
steadily approach, the true solution, terminating the iterative sequence when our 
results are sufficiently accurate for our purpose. 


Important considerations for iterative methods are therefore (i) to find a convergent 
method, and (ii) to find a method which converges reasonably quickly, i.e., one which 
has a reasonable rate of convergence. Different methods will be needed in different 
circumstances, and we shall rarely be able to estimate in advance the amount of 
computation required to achieve some specific accuracy. We shall, however, be'able 
to perform some useful analysis for the particular types of linear equations which 
arise from finite-diflerence methods applied to our elliptic boundary value problems. 


One main feature of such equations is the sparseness of the associated matrix, that is, 
the large number of zero elements in the matrix. This is manifest in the following 
matrix, typical of our problems, in which the nonzero elements are marked with a 
cross. 


x x x 
X x x x 
X x x x 
x R x m 
x x x (1) 
| 3 x X x x x 
MAN J I den) 


x X x 
x | X x x 
x X x 


For such matrices direct methods have the disadvantage that they use a large amount 
of computer storage. Any attempt to take advantage of the many zeros tends to bc 
frustrated by the fact that, in the process of elimination. few of the zeros within the 
dotted lines in the diagram can be preserved. They get "filled in" with nonzero ele- 
ments, giving rise to more arithmetic than would seem to be necessary. If interchanges 
are necessary to avoid induced instability then even more nonzero elements are 
introduced in other parts of the matrix. 


Iterative methods do not destroy the sparseness. and for the type of matrix shown 
above we can find some useful methods with satisfactory rates of convergence. In 
discussing these methods we shall assume that our linear equations are well-con- 
ditioned. They can be written in the matrix form 


Au = b, (2) 
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which is well-conditioned if small changes in the elements of A and b do not cause 
large changes in the solution vector u. This is clearly a reasonable assumption if 
our given problem is properly posed and our finite-difference equations are such that 
we get nearer and nearer to the true solution as the mesh spacings are steadily 
reduced. 


We know, moreover, from Unit M201 8, that the Gauss elimination method with 
pivoting does not suffer from induced instability, the solution obtained being the 
exact solution of a problem Au = b with small perturbations in the elements of A 
and b. If the problem is well-conditioned our computed result must therefore be close 
to the exact solution of Au = b. We can show, though we shall not prove this, that 
our iterative methods likewise do not suffer from induced instability. 


SAQ 5 


(a) Put crosses for the nonzero elements in the matrix (1) after two steps of Gauss 
elimination without interchanges. 


(b) Put crosses for the nonzero elements in the matrix (1) after two steps of Gauss 
elimination, where the "initial" matrix has rows | and 5 interchanged prior to 
the first step. 


(Solution on p. 38.) 


Notes 


(i) S: page 144, line 3 

A band matrix is one in which all nonzero entries are confined to the leading 
diagonal and those diagonals close to the leading one. That is, if a; ; is that 
element of a matrix A which lies in row i and column j, then A is a band matrix 
with band width 2k + 1 if 

aij = 0 for |i — j| > k. 
For example, a tridiagonal matrix is a band matrix with band width 3 since 

a; = O for |i — j| > 1. 
The matrix (1) exhibited earlier in this section has band width 9. 


(i) S: page 144, line — 2 to page 146, line 4 
Gauss elimination with pivoting is the subject of Unit M20/ 8. The idea of 
triangular decomposition of a matrix was also discussed in that unit 

(ii) S: page 146, lines 7 and 8 
One of the iterative methods we shall study later in this unit is the successive 
over-relaxation method (SOR for short) in which there is a parameter called the 
over-relaxation factor. The rate of convergence of the SOR method depends on 
this factor and as a result a large proportion of the theory associated with this 
method deals with finding a value of the parameter (called the optimum over- 
relaxation factor) which maximizes the rate of convergence of the method. 

(iv) S: page 146, lines 9 to 11 
A system of simultaneous equations can be written in matrix form as 


Ax —b 


where x = (x4, X5,....X,) is the column vector of unknowns. If we have an 
approximation X to x then the vector 


r=b-— Ax 


is called the residual vector of the given approximation. We can interpret the 
residual vector as a measure of how close the approximation x comes to satisfying 
the given system of equations. 


For a discussion of residuals and their uses see Section 8.4.1 of Unit M201 8. 
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11.3.2 Theory of Iterative Processes 


In this section we shall investigate iterative processes in some generality. This will 
enable us to find a general result about the convergence of iterative methods which we 
can then apply to the special methods used for the solution of boundary value 
problems. The ideas we shall begin with should be familiar to you from work done 
in both the Mathematics Foundation Course and the Linear Mathematics Course. 
In particular we adopt here the approach of Unit M100 28, Linear Algebra IV; 
you might like to review Section 28.2 before continuing with the present work. 


We begin by considering the system of equations written as 

Ax = b. (3) 
In general we can rearrange this equation as 

x = Gx + Hb. . (4) 
For example, adding x ~ Ax to each side of Equation (3) yields 

x=(I~ A)x +b (5) 


where 7 is the identity matrix. Equation (5) is a special case of Equation (4) where 
G =I — A and H = I. We shall return to Equation (5) later, but for the moment 
let us look at the iterative scheme which we can obtain from Equation (4) written as 

x! = Gx" 4 Hb, (6) 
Beginning with a certain initial approximation x, which can be chosen arbitrarily, 
Equation (6) yields successive approximations x (1 = 1,2,...). We would like to 
determine the.conditions under which we can be sure that the successive approxima- 
tions x converge and that the limit is a solution of the given Equation (3). The latter 
requirement is dealt with in our next theorem. 


THEOREM 3 

If the process defined by Equation (6) converges then lim x™ is a solution of Equation 
(3). 

Proof 

Suppose that 


lim x™ = 


Proceeding to the limit in Equation (6) we obtain 
x* = Gx* + Hb. 


That is, x* satisfies (4), and hence Equation (3), since (4) is simply a rearrangement 
of (3). 


It is clear that if we choose x* as the initial vector x! then all subsequent approxima- 
lions x (n = 1,2,...) will be equal to x*. This fact leads us to call x* a fixed point 
of the iteration. Whilst we are on the subject of terminology, we call iterative schemes 
of the type (6) stationary processes whenever the matrices G and H remain constant 
throughout the whole computation. ( Non-stationary processes have been constructed 
with the general form 


x™ = ginyti-1) 4 pun 


where the matrices G™ and H™ depend on n. In this course we shall deal only with 
stationary processes of the type (6).) 


Turning now to the conditions under which Equation (6) converges we state the 
following theorem. 


PDE! 


THEOREM 4 
For a stationary iterative process given by 


x= Gx") + Hb 


and sufficient that all eigenvalues 


jniti AHO SL , 
to converge for any initial vector x^ it is necessary d 
t | radius of G must be 


of the matrix G be less than one in modulus. That is, the spectra 
less than unity. (See Unit 8 for the definition of spectral radius.) 


We shall not prove this theorem for an arbitrary m x m matrix G. Instead. we shall 
verify it for the case when G has m linearly independent eigenvectors. The theorem is 
also true for a general square matrix. but the proof takes too long even though it 
requires no knowledge of linear algebra beyond that contained in M201. In addition 
we shall simplify matters by considering a process derived from 


Ax =b 


only for the case when A is nonsingular, so that the solution is uniquc. 


Proof 

Consider the process 
x? = Gx" 7" + Hb. 

and let x* be the solution of Equation (3), so that 
x* = Gx* + Hb. 

Subtracting the previous equation we obtain 


x* — x? = G(x* — xl), 


The vector e™ = x* — x" represents the error* in the nth approximation Lo 3 
and the last equation gives 


e = Gel!) = GTI Lo... = gne) 
If we suppose that {x™} converges as n ~ zc. its limit is x* by Theorem 3 and it 
follows that lim e = 0 as? ^x oc ; hence 

lim G"e'? = 0, 
Since this must hold for any initial vector x'®, and so for any error vector e, it is 
necessary that G" approaches 0. the zero matrix. Conversely if 

lim G" 


naa 
is the zero matrix it is clear that {x“! converges to x*. 
We now assume that the matrix G has m linearly independent eigenvectors which 


therefore may be used as a basis in the m-dimensional space of column vectors. 
Therefore, we may write 


where the v; are the eigenvectors of G and the c; are scalars. Since 
erm = G" 
we see that 


m 
e" — Y Gy. 
tel 


Now, by the definition of an eigenvalue we have 


my n oin 
Gy, = Avi 


* In Unit M100 28 the error vector was defined as the negative of this one. 
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where 4; is the eigenvalue corresponding to v;. Therefore, for lim e = 0 as n ~% 
we require 


m 
e" = J cv; approaches 0 as n ~ cc, 
i=l 


which can happen for arbitrary e% if and only if [4| < 1 for all i= 1,2,...,m. 
Hence, lim G" = 0 if and only if |2;| < 1 for all i, that is, if the spectral radius of G 
is less than unity. 


Tying up the two parts of the proof, we have proved our theorem for this special case. 
SAQ6 
Use the simple iterative scheme given by Equation (5) to solve the system of equations 
05x,— 08x, = 17, 
—0.2x, + 0.5x, = —0.5, 


with the initial approximation (x,, x2) = (0,0), performing the first five iterations 
only. What can you say [rom your results about the effectiveness of the simple iterative 
scheme in solving the given system of equations? Verify that the vector (5,1) is a 
fixed point of the iteration. 


(Solution on p. 39.) 


SAQ 7 
Can the system of equations — 
XQ x2 =4 
3x, + 7x, = 20 
be solved by the iterative method given by Equation (5)? 


(Solution on p. 40.) 


11.3.3 Convergence Rates 


We have seen in SAQ 6 that the simple iterative scheme 
x" = (J — Ax +b 


can converge very slowly. Since we can construct many different iterative schemes 
of the type 


x" 2 Gx"7! + Hb 


by choosing different matrices G and H it seems reasonable to choose the one which 
is most rapidly convergent for a particular problem. In order to be able to compare 
schemes we need a measure of the rate of convergence of any given choice. 


Once again if we assume that the matrix G has m linearly independent cigenvectors 
v; corresponding to eigenvalues 4; we obtain. for some choice of constants c;. the 
relationship 

m 

e" — LY cv; (7) 

i=) 
where e"" is the difference between the true solution of the system and the nth approxi- 
mation as defined in Section 11.3.2. For simplicity we shall assume that 4, is the 
unique eigenvalue of largest modulus, ie. || > [44 for all i # 1: thus the ratios 
(4,/4,)" approach zero as n becomes large. Therefore, writing (7) in the form 


; 43^ A" AJ 
e" = 2 40v, + [2] esa [I eva o e E] eun 
Ay Ln ^l 
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we see that 
e? ~ Mey for n large. 
Similarly 
et Dw Bley, for n large. 


and therefore e"*! ~ A,¢ for sufficiently large n. As |A,|, by definition, is the 
spectral radius p of G, we can say that the error in the approximation ultimately 
decreases approximately by a factor 1/p in each subsequent iteration. 


Suppose that the error in the kth component of x" is e. which is given by the approxi- 
mation 
£ = lee = ical = p"lel say, 


where n, , is the kth component of the eigenvector v, associated with the eigenvalue 
A, and $ = cyX,.,. Taking logarithms we have 


_ leg (ple), 
-logp 


Hence n, the number of iterations required to reduce to c the error in each component 
of the solution vector, is inversely proportional to —log p. We therefore define the 
asymptotic rate of convergence of the iterative method as —log p, and in the choice 
of a good iterative method our aim is to select a matrix G whose spectral radius is 
as small as possible so that the rate of convergence is as large as possible. 


SAQ 8 
What is the asymptotic rate of convergence of the scheme used in SAQ 6? 
(Solution on p. 40.) 


We have assumed in our analysis that |2,| > |2;| for i + 1. Since G is a real matrix 
its complex eigenvalues occur in conjugate pairs, which must have the same modulus; 
hence 4, is real. In this case we can approach the question of the rate of convergence 
somewhat differently and in the process compute an approximation to 2,. Suppose 
that A” is the correction, or displacement vector, to x™ at the (n + 1)th iferation. 
That is, we define 


Am = yt yon 
The iterative scheme can then be regarded as the summation of the infinite series 
* ym p AD) AED Ee, 


Now we have seen that for sufficiently large n the error vector e" = x* — y 
satisfies the relation 


etna A, e"! 


and it is easy to see, by similar reasoning, that the displacement vector satisfies a 
similar equation given by 


AMT Y x AC 
(This can hold only if 2, is real.) Thus, for sufficiently large n 
x*-x x"! + AML HA, + AP + iie) 
and if [2,| < 1 we can sum the infinite geometric series and obtain 
AM 
IA (8) 


This result tells us that the an NU 
URL current correction A? should really be multiplied by 


x* c xU 4 
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This is an important result because if we want an approximation to x* with error no 
greater than c we might be tempted to terminate the iteration at the first n for which 
Al, < £. where the || | notation denotes the uniform norm, 


Al] = max Aj. 


If 2, is 0.99, for example, which is not uncommon, such a termination would give a 
very poor result, and the iteration should be continued until A = 0.01 c. In 
general we do not know 2, but it can be estimated from the ultimate ratio of the 
corresponding components of successive displacement vectors which we can easily 
compute. Such information is very important because it not only tells us when to 
stop the iteration but it can also help us to accelerate the convergence. This is demon- 
strated in the following SAQ. 


SAQ9 


From the results of your iteration in SAQ 6, estimate the eigenvalue of largest 
modulus of the iteration matrix. Using this estimate, try to compute a better approxi- 
mation to x* from the approximations x'^ and x? : 


HINT: Consider Equation (8). 


(Solution on p. 40.) 
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11.4 ITERATIVE METHODS FOR BOUNDARY VALUE 
PROBLEMS 


11.4.1 Jacobi and Gauss-Seidel Methods 


This section is intended as an introduction to one class of iterative schemes arising 
naturally in the solution of the finite-difference equations which approximate elliptic 
partial differential equations. The methods we are describing here can also be used 
for the solution of initial value problems. For these, however, direet methods of 
solution (Gauss elimination) are far more useful, and it is for boundary value problems 
that iterative methods become really important. 


Throughout Section 11.4 we shall investigate the solution of Poisson's equation 
QU PU 
Il dut. 
Ox? oy? . 
in a rectangular region. We have rewritten Poisson's equation in this form to conform 
with the usage in the next reading passage from S. Furthermore we shall discretize 
the problem using a square mesh with equal mesh spacings h in both the x- and 


y-directions* We shall assume that Dirichlet boundary conditions are specified, with 
U given at all points on the boundary. 


+f=0 


In the following reading passage many references are made to passages in Chapter 3 
of S which you have not read. You need not refer to Chapter 3 at all since the notes 
in this text contain all the extra information you will require. 


Notes 


(i) S: page 147. Jacobi method 
If we write the Jacobi method in matrix form we obtain 


wD = Ju 4p 


where J is the square matrix of order (p — 1)(4 — 1) given by 


[BI 
I B I 
I BI 


NS 
ll 
He 


in which each block is a square matrix of order (q — 1), with 


* This presumes that the ratio of the lengths of adjacent sides of the domain is rational. 
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and / the identity matrix. The column vector b is determined by the /; j and b, ;. 
We call J the Jacobi iteration matrix. We can see immediately that the iterative 
scheme is in the form of Equation (6) of Section 11.3.2 and therefore to investigate 
its convergence we need to find the spectral radius of the matrix J. 


In fact, as shown in the Appendix, the cigenvalues 4 of J are given by 


dij 


l in jn : . 
E + 2 cos!) f=1,2,....p-yj=1,2.....¢-1, 
and the spectral radius p is therefore 
i 
pV) == [ss + es? . 
2 p q 


Hence p(J) < 1 (since p,q > 1) and the Jacobi method is convergent. 


S: pages 147 to 149, Gauss-Seidel method 
In matrix form Equation (5.7) of S: page 148 becomes 


wt) 2 [y^*D 4. Uu? + p, 


where 
B, ro = 
L B; 10 
1 B, 
L-i > B= " 
I B, 10 
L I B,] L 1.04. 
and 
M ] Fo 1 J 
"er 0 1 
By I 
U-i » B= 
By I ee 
L By J d oJ 


Land U are strictly lower and strictly upper triangular matrices respectively of 
order (p — 1)(q — 1), B, and By are square matrices of order q — 1 and J is the 
unit matrix of order q — 1. Notice that L + U = J, the Jacobi iteration matrix. 
The Gauss-Seidel method for our problem can be rewritten as 


(I — Lu"*? = Uu" + b 
and we see that 
wet) = (I — L)^! Uu'? + (I — L)!b 


which is of the general form of the iterative schemes dealt with in Section 11.3. 
We call the matrix (Z — L)^!U the Gauss-Seidel iteration matrix. Fortunately. 
for the matrices which we are considering, it is not necessary to find the eigen- 
values of the Gauss-Seidel iteration matrix directly since it is possible to express 
them in terms of the eigenvalues of the Jacobi iteration matrix. We shall not do 
this here. However, notice the relationship, given at the top of S: page 149. 
between the spectral radii of the two iteration matrices. 
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General Comment 

We may use an iterative scheme to solve the matrix equation 
Ax =b 

where 
A=I-L-JU, 


with J the identity matrix and L, U being strictly lower and strictly upper triangular 
matrices respectively. In general we call the scheme 


xt Py + Ux" +b 
a Jacobi method, and 
xO n pxe*D p Ux™ 4b 


a Gauss-Seidel method. 


SAQ 10 


The five-point formula, when applied to Laplace's equation for the region shown with 
Dirichlet boundary conditions, yields the system 


Wr dus ~ =b’ 
aly + — dus = by, 
— dur + u = by. 
(a) Determine the Jacobi iteration matrix J. 
(b) Determine the Gauss-Seidel iteration matrix G. 
(c) Verify that the spectral radii are related by 
[oN]? = p(G), 
for this problem. 


(Solution on p.40.) 


SAQ 11 


(a) Show that the Jacobi method is convergent and the Gauss-Seidel method diverg- 
ent when applied to the system - 


uU, + Quy — 2u; = 1, 
Wy + ty + 3 = 3, 
2u, + 2u, + us = 5, 


. PDE 11.4.1/11.4.2 
(b) What do you find when you perform the Jacobi iteration with a starting value 
given by 
a 
x =| B 
y 
and exact arithmetic is used? How could you have predicted your result? 


(Solution on p. 42.) 


11.4.2 Accelerating the Rates of Convergence 


Consider the finite-difference solution of Poisson's equation in the region shown in 
the figure when Dirichlet conditions are given on the boundary. 


If the mesh lengths are equal, the five-point formula produces a set of linear equations 
of the general form 


uy = Hina + Ui 1j t tij- + tijt) + bij 
and typical equations relevant to the figure are 
uy = dius + ug) + by, 
ua = A + Uy + us) + ba, (1) 
Us = (us + Uy + Ug + ug) + Ds, 
etc. 


Now the essence of an iterative method is that we scan the mesh points systematically. 
Having reached point r we obtain a new value for u, using the relevant equation from 
the set (1). With this new value of u,. some of the other equations fail to be satisfied. 
whence the need for iteration. 


The Jacobi and Gauss-Seidel methods perform the operations slightly differently 
from each other. Suppose that, for example, we start with some initial guessed values 


up um. en ug’ relevant to the figure. At point 1 both methods produce the new 


approximation u!! from the formula 
ul? = iut + ul) + by. 


However, at the next stage, the adjustment at point 2, the Jacobi method uses the 
formula 


ut Hut? + lY + ul) + by, (2) 
while the Gauss-Seidel formula is 


ul! z ia! +? + nl) + ba. 


Thus. for subsequent points. we do not use the new (and hopefully better) valuc up 
in the Jacobi method. whereas in the Gauss-Seidel case we take advantage of this 
better available value in order, again hopefully, to accelerate the convergence of the 
process. If we perform the adjustments scanning the points in the natural nier 
ordering shown in the figure we see that, for example. in the Gauss-Seidel metho: 
the adjustment at point 5 is given by 


3 
ul) = lu + a) + uO + uP) + bs, (3) 
which uses the two “better” values already known for points 2 and 4. 


The rate of convergence of the Jacobi iteration is clearly independent of the order in 
which the mesh points are scanned since the arithmetic is independent of the ordering. 
For the Gauss-Seidel method, however, not all orderings will involve the same 
arithmetic, For example, suppose we scan the points in the figure in the order 


1,2, 3; 6, 5, 4; 7, 8, 9: 
then the adjustment at point 5 is given by 
uf? = Md! + ul + uf? + uP) + bs e 


since we already have the new values at points 1, 2, 3 and 6 but not at 4, 7, 8 and 
9. Equation (4) is obviously different from Equation (3) and we might expect that one 
of the orderings will have a faster convergence rate than the other. 


In discussing this topic it is convenient to consider at the same time the idea of 
introducing some parameter into the matrix equations, hoping to obtain even faster 
convergence with the choice of a suitable numerical value for this parameter. 


The motivation for this was the discovery by early workers in this field that quite 
commonly all the successive changes in any pivotal value in successive steps of, say, 
the Gauss-Seidel process, are of the same sign. For example, approximations at a 
particular point may have successive values like 100, 120, 130, 135, 137, 138, converg- 
ing monotonically to the solution. (SAQ 6 is a typical example.) This indeed must 
happen eventually if the eigenvalue 4, of largest modulus in the iteration matrix is 
real and positive since, as we saw in Section 11.3.3, successive changes are ultimately 
connected by the equation 


NUES 24,A0 


for sufficiently large n; in practice, a positive 4, turns out to be the rule rather than the 
exception. In such cases one might accelerate the convergence by making a larger 
change A™ than is needed to satisfy (temporarily) the relevant equation of the set (i ). 
Now, if the original equations are given by 


Au =b 


where A = I — L — U, L being strictly lower triangular and U Strictly upper tri- 
angular, then the Gauss-Seidel method is given by 


wit) = Lott + Uu? + b 
and the current correction is just 
AU x nt D og) ope pulito (U = Du, 


If we multiply this correction by c, where usually 1 < c 


< 2, the iterative scheme 
becomes 7 


wet ee un = «ib + Lut) + (U - Du 
or 

U — oLy"*!" = i wl «Uu? + wb. (5) 
This is the iteration equation for the successive over-relaxation or SOR method 


("over" because at each stage we are "overdoing" things). Notice that the SO 
. a 
method reduces to the Gauss-Seidel method when w = 1, 5 
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The two questions, of choosing a suitable ordering of the equations and a suitable 
value for the parameter w, have received much attention in the last few years. and 
in the next section we give an introduction. without detailed proofs, to the more 
important results of this research. 


11.4.3 Consistent Ordering, Property 4 and SOR 


We consider further Poisson's equation in the square domain of Section 11.4.2. 


Suppose we order the points for the iterations of Gauss-Seidel or SOR methods 
along the diagonals, that is, in the order 


1:4,2:7, 5, 3:8. 6:9. 
In matrix form the equations become 
ut» = Lut) + Uu? T b, (6) 
where 7 — L — U has the form 
|] 4 2 7 5 3 8 6 9 


x |x x 
x | x x x ~ 
x x x x 
j ct 
x x x 
x x x x x 
x x x 
Ix x x x 
X x Xx 
x xix] 


in which the crosses represent nonzero entries. Equation (6) implies that we sc: 
mesh in the diagonal order quoted, adjusting successive values from equations like 
(1) of the previous section, but using the latest available pivotal values on the right- 
hand side of (6). For example. Equation (3) for the Gauss-Seidel method becomes 


ue! n o duet A uP O ug + ug bs. (7) 


which is of the same form as the natural ordering 1.2. 3..... 9 because in cach case 
uy and u4. but not uy and ug. have already been "adjusted". 


Consider next the ordering 


1.3.5. 7,912. 4. 6. 8. 


: ; atex Te ' has 
the so-called black-white or checker board ordering. The matrix 7 — L — U now ha 
the form 


x x x 
x x x 
x [-* x x x 
x x x 
x x x 
x x x x 
x x x x 
x x x x 
L x x x x J 


and in this case Equation (7) becomes 
Wg) d ae ut ae ub? + up) + bs 


since none of the pivotal values on the right-hand side has yet been adjusted. When 
We come to compute a new u$! on the other hand, all the other pivotal values in the 
relevant equation have the superscript (n + 1). 


The partitioned matrix of Equation (6) with the diagonal ordering has the diagonally 
block-tridiagonal form 


[D. F, 
E, D, F 
s (8) 
En-2 Dm- Fray 
where the D; (i = 1,2,...) are diagonal matrices. not necessarily of the same order. 
The partitioned matrix obtained from the checker board ordering. which we can 
write as 
[D, F 
LE Dy. ` (9) 


is of special diagonally block-tridiagonal form. 


SAQ 12 


What is the form of the matrix obtained using the natural ordering 1,2,3; 4,5,6: 
7,8,9? Is it of either of the special forms (8) or (9)? i 


(Solution on p. 43.) 


It turns out that there are other orderings, called consistent orderings, which do not 

produce matrices of types (8) or (9), but for which the arithmetic of the Gauss-Seidel 

or SOR method is identical with one of these forms. That is, at each stage of the 
gl 


iteration the values obtained at each point are identical with those obtained using 
one of these forms. P 


For the diagonal ordering the following table shows, for the computation of succes- 
sive u{"*) those values in the respective equations which already h 
(n + 1) and those which, not yet adjusted in the current iterati 
superscript (n). 


ave the superscript 
on cycle, have the 
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points with superscript (n + 1) points with superscript (n) 


ral none 4.2 
4 1 7,5 
2 I 3,3 
7 4 8 
5 4.2 8.6 
3 2 6 
8 45 9 
6 5.3 9 
9 8.6 none 


Here the numbers in the second column refer to the elements below the diagonal of 
the matrix corresponding to the row number given in the first column, and the third 
column refers to the corresponding elements above the diagonal of the matrix. The 
next SAQ shows, among other things, that there is another ordering not in the form 
of (8) but consistent with that of the table above. 


SAQ 13 


(a) Why are the statements about orderings for the SOR method equivalent to 
those made about the Gauss-Seidel method? 


(b) Show that the natural ordering of SAQ 12 is consistent with that of Equation (6). 
(c) Show that the ordering 
1, 2, 3; 6, 5, 4; 7, 8, 9 
is not consistent with either the diagonal ordering or checker board ordering. 
(Solution on p. 43.) 
We now give two definitions. 


A square matrix is said to have Property 4 if, by transposing pairs of rows and cor- 
responding columns, it can be transformed to the form of the matrix (9). In particular, 
the matrix (8) has Property A. (To see this, rearrange the matrix so that the order of 
diagonal blocks is Dj, D3, Ds,..., Dz, Dy...) 


Let Q be a matrix obtained from P by successively interchanging pairs of rows and 
corresponding columns. For example, let Q be obtained from P by interchanging 
the ith and jth rows, and interchanging the ith and jth columns. Then, if Ej; is obtained 
from the identity matrix by interchanging the ith and jth rows, 


Q = Ej "PE. 


In general, if T is the matrix obtained from the n x n identity matrix by reordering 
the rows to 


TEETETTIIN 
then 

Q =T PT 
is the matrix obtained from P by reordering its rows and columns to the new order. 
It is clear that the main diagonal of Q consists of the elements of the main diagonal 
of P in the new order. We say that Q is consistent with P if the set of nonzero elements 


below (above) the main diagonal of Q is a permutation of the nonzero elements 
below (above) the main diagonal of P. 


In the context of the Gauss-Seidel method the matrix obtained after reordering 
I — L — U is consistent with I — L — U if the orderings are consistent in the sense 
discussed previously. If] — L — U has Property A then, for brevity, we shall say that 
I —L — U is consistent if it is consistent with any matrix of the form (8) or (9). 
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The point of making our definitions is that if 7 — L — U has Property A, the asymp- 
totic rate of convergence of the Gauss-Seidel or SOR methods is the same for all 
orderings such that the matrix J — L — U is consistent. It is, of course, not surprising 
that, for example. the diagonal ordering and the natural ordering should yield the 
same rate of convergence, since we showed in SAQ 13 that these orderings are 
consistent with each other and give rise to precisely the same arithmetic. The more 
remarkable thing is that the diagonal ordering and the checker board ordering, which 
are not consistent with each other and give rise to different arithmetic, should still 
have the same asymptotic rates of convergence. It is the Property A quality of the 
original matrix which produces this result. 


It might be asked why we concentrate on consistent orderings and on matrices with 
Property A. First, it can be shown that the asymptotic rate of convergence for all 
inconsistent orderings is slower than that for consistent orderings. Second, for 
matrices with Property A. and for all consistent orderings, we can find an optimum 
value for the parameter which gives the fastest rate of convergence for the SOR 
variation of the Gauss-Seidel method. 


For this purpose we need to investigate the spectral radius of the SOR iteration matrix 
Sa = U- oL)! — wI + wU} 
obtained from Equation (5). 


First we shall prove a lemma. (Omit the proof if you are short of time.) 


LEMMA 


If the matrix J — L — U, where J is the identity matrix and L(U) is strictly lower 
(upper) triangular, has Property A and is consistent then the eigenvalues of 
pL 4- 7! U (u + 0) are independent of p. 


Proof 
Since J — É — U has Property A and is consistent, there exists a matrix T'such that 


Hh Fy 


E h Fi 
T"(-L-UT- 


En-2 du. F, 


m-i 


E I 


where the elements of E,, E;...., £,,-, are the elements of L, and the elements of 
F,F. F,-, are the elements of U, suitably reordered. (The diagonal blocks 
are all identity matrices. Why?) Now, it is clear that multiplying J, Lor U by a constant 
will reproduce this multiplication in precisely those elements of the block matrix 
which come from J, L or U respectively. Hence i 


m-1 "m 


Al, a! F, 


TAL — uL — u^! UT = 


Em- AL, 


Now the eigenvalues of uL + 7 'U are the roots of the equation 


MI — uL- po 'Ul = 0. 
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Since a similarity transformation preserves the eigenvalues we need look only at the 
matrix on the right-hand side. Suppose J; is of order pj. We leave the first pi TOWS 
and columns as they stand. Then we multiply the next p; rows by j^! and the next 
p; columns by 4, obtaining the matrix 


Al, F, 


WE, Ay g^! Fs 


HEm-1 Alm 


Proceeding in this way, we multiply the next P3 rows by 17? and p, columns by 2, 
etc, finally multiplying the last p, rows by 1^ "*! and the last Pm columns by y^^! 
leaving us with the result that ^ 


A F, 


ll 


A — uL — u^! u| 
En-2 Aly-i Fm- 
Ej. AL, 


= |T — L — U)T| 
= |Al — L— UI, 


since the row and column operations have, in total, multiplied the determinant by 1. 
Hence the eigenvalues of yL + 47 'U are the same as those of L + U; ie, they are 
independent of y. 


We now use this result to relate the eigenvalues of the SOR iteration matrix to that 
of the Jacobi iteration matrix. 


If 4 is an eigenvalue of S, then 
S, — All 2 0. 

Then 
(2 — oL) U + w(U — 1) - al] =0 
so that, on multiplying by |/ — wL], we obtain 
Ic-o(U-1-20-oLIz0 
since determinantssatisfy the rule 

PQ| = |PIIQl. 

Hence 


|XU + 4L) - G + w — 1| = 0, 


ie. 
A3oGL + 2730) - (A+  — 1| = 0. 
Finally 
i zu 
[ab aes) Eug. (10) 
A W 


PDE 


, iste e eigen- 
Now by the Lemma. if J — L — U has Property A and is ARE Riu asa 
values of the matrix uL + u^! U are independent of p. It follows that t dedos 
of 33 L + 47 3U are identical with those of L + U, and these are precisely kable 
values « of the Jacobi iteration matrix. Equation (10) then gives the remarka 


relation 


Ac-to-l1 (11) 
Oe 
between the eigenvalues a of the Jacobi iteration matrix and the eigenvalues 4 of i 
SOR iteration matrix. With w = 1, incidentally, we confirm the result 2 = « whic 
we have already observed, for example, in SAQ 10 and which guarantees, when 
la] < 1 for the convergence of the Jacobi method, that the Gauss-Seidel method 
converges twice as fast. 


Proceeding from Equation (11) we can show that the value of the parameter a which 
minimizes the spectral radius of the SOR iteration matrix, thereby giving the fastest 
possible rate of convergence, is given by 


2 


SE esi 


where p(J) is the spectral radius of the Jacobi iteration matrix. Moreover, the spectral 
radius of the SOR iteration matrix is then given simply by 


P(54,,) = Wop — 1. 


We have developed the program $SOR 321 to solve up to ten equations iteratively 
by SOR; a listing follows. 


18 PRINT “THIS PROGRAM USES THE METHOD OF SUCCESSIVE OVER RELAXATION 
28 PRINT "TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO SIX DEC.PL." 
38 PRINT 

48 DIM ACIO,112»XC19J,ZC102,A8C31 

58 PRINT "TYPE NUMBER OF EQUATIONS(«-10)''; 

60 INPUT N 

76 PRINT 

8Ø PRINT “TYPE ELEMENTS OF COEFFICIENT MATRIX»ROW BY ROW" 

98 FOR I=1 TON 

180 PRINT "ROW'";Ij 

118 MAT INPUT ZCN) 

120 FOR J=1 TON 

130 LET ACIsJ1=Z0d] 

140 NEXT J 

150 NEXT I 

166 PRINT 

178 PRINT "TYPE R.H.S. CONSTANTS" 

180 MAT INPUT ZENJ 

198 FOR I=1 TON 

288 LET ACL,N*11eZLIJ 

2180 NEXT I 

220 PRINT 

238 PRINT "TYPE INITIAL APP " 

248 MAT INPUT ZEN) ree aay 
250 PRINT 

260 PRINT "TYPE OVER RELAXATIO 
278 INPUT w 

280 MAT X-ZERINJ 

290 PRINT 

309 PRINT 

310 LET C=O 

320 FOR I=1 TON 


N FACTOR"; 
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338 LET M-ACI,IJ 

348 FOR J=1 TO N+1 

350 LET ACI»>JI=ALI,JI/M 
368 NEXT J 

378 NEXT I 


388 LET C=C+1 

398 LET E=@ 

400 MAT X-Z 

416 FOR I=1 TON 

426 LET S1=S2=ø 

436 FOR J-I«1 TON 

440 LET SI-SI*ACISJIX*XCJJ 


450 NEXT J 

460 FOR J=1 TO I-1 

470 LET S2=S2+ACIsJI*ZCJ) 

488 NEXT J 

490 LET ZEIJeCI-WO*XCIIHW*CACIoN*12J-S1-S2) 
500 LET K-ABSCZCIJ-XCIJ) 

518 IF K <= 5.E-07 THEN 539 

528 LET E=1 

538 NEXT I 

540 IF E=} THEN 389 


550 PRINT 
5680 PRINT "SOLUTION" 

5780 PRINT "-------- 

588 MAT PRINT Z 

590 PRINT 

696 PRINT "NUMBER OF ITERATIONS NEEDED IS'';C 
610 PRINT 


628 PRINT "REPEAT WITH ANOTHER OVER-RELAXATION FACTOR/INITIAL APPROX. 
630 PRINT "YES OR NO''; 

640  INPUT A$ 

650 IF ASz'"YES" THEN 239 

660 END 


SAQ 14 
Consider the equations of SAQ 10, and its solution which gives 


[p = $ = p(G). 
(a) Find the optimum w for the SOR method and the spectral radius of the SOR 
iteration matrix. Does the original matrix have Property A? 


(b) Determine the improvement in convergence by performing both Gauss-Seidel 
and SOR iterations with the optimum value of w, for the problem of SAQ 10 
with b, = L5, b; = 3, b = 3.5. You may use the library program $SOR321. 


(Solution on p. 44.) 


Two questions remain. Firstly, to get w,,, we apparently need an estimate of p(J). 
This is sometimes known theoretically, as in note (i) of Section 11.4.1. Alternatively 
we iterate for some time with the Gauss-Seidel method, computing the spectral radius 
of its iteration matrix, which is [p(J)]^, by the method suggested at the end of Section 
11.3.3. In practice, it is better to overestimate the spectral radius of the Jacobi iteration 
matrix because it will yield a smaller decrease in the rate of convergence of the SOR 
method than would a comparable underestimate. 


Secondly, what can we do if our matrix does not have Property A? In such cases the 
SOR method has been used with some success, though without theoretical justifica- 
tion, and much current research is devoted to this problem. 
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PDE 11.5 


11.5 SUMMARY f 


We began the unit by investigating two simple applications of the finite-difference 
method applied to boundary value problems. We saw that the numerical method 
requires investigation in two different areas. The finite-difference method gives rise 
to (large) systems of linear algebraic equations and we need to know the conditions 
under which the solution of the equations converges to the true solution of the 
boundary value problem. We were able to prove, for bounded regions, that the 
five-point formula for the Laplacian is convergent and moreover we were able to 
obtain a bound on the error of our approximate solution. Both of these results are 
consequences of a maximum principle which we also established. 


The second half of the unit concentrated on the problem of solving the systems of 
algebraic equations which arise in the finite-difference method. In particular, we 
illustrated that iterative methods have advantages over other methods for our 
particular problems. We investigated iterative methods and showed that in general 
the iteration matrix must have a spectral radius less than unity, and that the smaller 
the spectral radius the faster is the asymptotic rate of convergence of a method. 


In the final section we looked at some particular iterative methods for boundary 
value problems, namely the Jacobi, Gauss-Seidel and SOR methods. We found that 
the SOR method could be optimized to give the fastest rate of convergence of these 
methods. However, our analysis depended upon the initial matrix having Property A 
and the mesh points being ordered in a manner consistent with this property. Under 
these circumstances we were able to state a result giving the optimum relaxation factor 
yielding the fastest rate of convergence. 


An alternative treatment of some of the material presented in Sections 11.3 and 11.4 
may be found in S: pages 76 to 79. 
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PDEII SAQ 1-2 


11.6 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 


We are given 


a a 

o 357 = -16 -1<Sx<1-1<y<1, 

ull, y) =0 -1€y&l, 

ài 

Fy Oe = cuis 1) -1€x<1, 

Putting y = 3+, we obtain 

ĝu au 1 = 1 

ax? t ge ~16 ERILILICIPII-. 
" 1 

u(1, 5) 2 0 -aS 

ðu 

a; 39 = —3tu(x, 31) -1€x€1. 


If we take a mesh spacing of 1 in the x-direction and aya in the y-direction, the central- 
difference formula is the same as that given on the last line of S: page 134, and the 
numerical results remain the same. 


Solution to SAQ 2 


See the solution given on $: page 162, lines 5 to 8. The following figure may be of 
assistance. 
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PDE ll SAQ 3 


Solution to SAQ 3 


The differential equation is approximated by 
1 
gp T)0)u,- 324,20 


and the derivative boundary condition along x — 1 by 


1 1 

ALIS = uis.) =F 2 Ui je 

The region is shown in the following diagram which also shows the mesh point 
numbering and the extension of the mesh to enable the derivative boundary condition 
to be dealt with by the central-difference formula. 


Zi 


ou _ Qu a 
aci rri 
> 
x 


We see that there is symmetry about the y-axis and therefore taking h = 1 there are 
35 unknowns to be evaluated. 


Typical equations are: 
for point 1, 
2u; + Ug — 6u, = 0, 
for point 2, 
uy + u3 + uy ~ 6u, = 0, 
for point 3, 
Uy + Ug d ug — 6u; = 0, 
for point 4, 
ug + Us + ug — 6u, = 0, 
for point 5 (two equations since it is a boundary point), 
Ug +u-4 cuo — bus = 0 
Unig — Ug = —Àus 
and eliminating u- 4 gives 


2u4 + Uio — Gus = 0. 
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Repeating this process for the remaining points yields 


five positions due to the boundary condition u — 1 along y = —1. 
Solution to SAQ 4 
The local truncation error T; of the five-point formula is given by 


_ Uii; Z 2U;j + Ui; $ Ui jad = 2U,;t U,; 


HTUE —: íj-l —]:; 
T, (Ax)? (Ay? "n 
By Taylor's Theorem 
0*U,; | (Ax)* 0*U,; 
Urs, 7 Uj + Uii m A n pL SE às t OAS), 
#U,; | (Ay)* 0*U,; 
Ui ja 7 2U + Ui j-i (y+ 12 à + O((Ayf). 
Hence, : 
9U,,  0^U,; 1 FLU , 9*U,; 
T XS D T fu 1 (Ax) Ae * (Ay) E 


+ O((Ax)*) + O((Ay)*). 
Since the differential equation gives 


au PU 
Bet BFF 
we find that 
1 2 9*U; 5d " (Ay 
- c L + (A + O((Ax)*) + O((Ay)*). 
T; i| a + A») 3 ((Ax)*) + 0((Ay)*). 


ooooocooo oo 


[-6 2 1 |] fy 

1-6 1 1 uz 

1-6 1 1 uz 

1-6 1 1 [" 

2 -64 1 Us 

I -6 2 1 Ug 2 

1 1-6 | 1 uy 

1 1-6 1 1 tig 

1 1-6 | 1 Ug 

1 2 -64 1 uzo 

L JL 
which gives the required form. Note that the constant vector b has —1’s in the last 
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PDE 11 SAQ 4-5 


Taylor’s Theorem with remainder gives 

ðU; ?QU,, PU he atU 
ig +h Uj ts CRM s+ 7 
Y Ox 2 ax? 6 óx 24 dx 


for some 0, e (0, 1), 


Urey = U (x; + Oh yj 


where h = Ax. Also, 
aU; Rh QU, WAU; hU 
Uu = Ug fU xa e B "xx 


for some 6, e (0, 1), 


(x; — 03h, yj) 


and so, 


2U,, h* aU 
Mb (x; + 03h, yj) 


Uji4,;—- 2Ui; + Ui ij ox IET 


for some 0, e (— 1, 1), 
using the Intermediate Value Theorem given in note (ii) of Section 8.2.1 of Unit 8. 
Similarly, on putting k — Ay, we obtain 


OU, . ktatu 


Uijay — 2U + Ui, = k? a + pg OO + 0,k) 


for some 0, e(— 1, 1). 
Therefore, 


QU, 2U; 
Ty = A * dy = f 


1 atu atu 
+ lense + Osh, y) + eae Y+ x 


The differential equation gives 


ou PU, 
ox) ^0y — 
and so 
1 atu o^U 
Thy -3[" Se + Ozh, yj) + eS oen + s]; 
Thus, 
1 atu au ` 
IT; | < al" 330 + Osh, y)| + k? aon + 04k) | 


1 
x —(h 2 
pM, + My) 


where M, and M, are the largest values of |@*U/dx*| and |@*U/dy*| in D U C, where 
C is the boundary of D. - 


Solution to SAQ 5 
In the following matrices the entries marked x are the nonzero entries. The entries 


marked ® are positions of nonzero entries where a zero entry occurred in the initial 
matrix. 
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(a) initial matrix 


K (b) initial matrix 
x x x ] x X x x 
x x x x x x x x 
x x x x x x x 
x x x x 
x x x x x x x 
x X x x x x x X x 
x X x x x x x x 
x x x x x 
^ x x x x x 
B x x x x x x 
x x x x x 
x x x 
- J L 
first step first step 
[x x x 7 [x x x x 
x x @ x x x @ x e 
x x x x x x x 
x x x x 
[-] x x x x x 9 e 
x x x x x x x x x 
x x xe x x x x x 
x x x x x 
x x x x x 
x x x x x x 
x x x x x 
x x x 
second step second step 
[* x x ] [ x x x x 
x x @ x x x @ x | 
x x|lo @ x x 6 ® x [:] 
x x x x 
[:] x x x [:] x @ @ 
@ x x x x ®@ x x x @ 
x x x x x x x x 
x x x x x 
x x x x x 
x x x x i x x 
x x x x x 
x x x J L 
Solution to SAQ 6 
For the simple iterative scheme given by 
x" = (I — Ax" ^" +b 
we have 
05 08 1.7 
POA) ga os SYS] os 
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Starting with x = (0, 0) we obtain 
afos ceno] T £7]. ui 
* T7lo» osjlo| [-05|" | -05]' 
m.[95 os] 17] f r7]. 2x] 
* 7lo2 os]|-os|*[-05; ^ | -041 
Similar computations give 
x?! = (2447, — 0.275), 
x = (2.7035, — 0.1481), 


x'5) = (2.93327, — 0.03335). 


The process should converge to the solution of the linear equations, x* — (5, 1) 
(unique, since det A + 0) since the eigenvalues of J — A satisfy 


(0.5 — A? = 0.16, 


so that 4, = 0.9, 2, = 0.1. But clearly the convergence is very slow. To verify that 
(5, 1) is a fixed point of the iteration we check that 


05 O08][5 n s] P C] 

02 05/1 -05! Ll 
Solution to SAQ 7 
For this simple iterative scheme, the iteration matrix I — A is 

0-1 

-3 -6 
whose eigenvalues 4 are given by A= —3 + J 12. The spectral radius of the iteration 
matrix is therefore greater than 1 and the scheme will not converge. 


Solution to SAQ 8 


The iteration matrix of SAQ 6 has eigenvalues given by 2, = 0.9 and A, = 0.1 (see 
solution to SAQ 6). Hence the spectral radius p = 0.9 and the asymptotic rate of 
convergence of the method is therefore — log 0.9 ~ 0.10. 
Solution to SAQ 9 
From the results of SAQ 6 we find 

A = (1.7, -0.5) — A? = (0.297, 0.135), 

A" = (0.45, 0.09), AG) = (0.2565, 0.1269), 

AM = (0.22977, 0.11475), 


The ratios of the components of A?! to A) are 0.864 and 0.940 and the ratios of the 
components of A to A?! are 0.896 and 0.904. We deduce that PES 
a better approximation than x? to be 


AO 
1-4, 
= (5.0012, 0.9994), 


which is very close to the true solution! 


0.9 and expect 


x? + =x™+ 104% 


Solution to SAQ 10 


The system can be written as 
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Hence L— 


eo 
e 


* o o 
o 
E 
a5 
a. 
iod 
LI 
o 
oO o s- 


oO Ae 
C Ae 


0 
J=L+U= 


Oo ae 
B 
oO ae 


(b) We have 


(c) 


I-L=|] —- 


o r 
Í 

a 
- 


so that 


(q-LyY!'- 


al 
E 


and 


G-(-L-'U-|0 + 


zb e 


0 gd 
For the Jacobi method, p(J) = max|a| where 


Hence 


and 


Therefore, 


ply) =. 


2/2 


For the Gauss-Seidel method, p(G) = max|2| where 


-À 1 0 
0 -4 i [|20 
0 d w- À 
Hence 
-44-4 =0 
and 
4=0,0,§ 


SAQ 10 


4l 


Therefore, 
p(G).= à. 
Clearly [p(J)]? = p(G). 


Solution to SAQII 
(a) The given system may be expressed as 
1 2-2||u, 1 


2 2 Vj} uy 
The Jacobi iteration matrix is 
0-2 2 
J={-1 0-1 
-2-2 0 
and its eigenvalues satisfy 
-4-2 2 
-1 —-À -1/=0 
-2 -2 -2 


which gives 4° = 0 and all the eigenvalues are zero. The iteration method con- 
verges since p(J) — 0. 


The Gauss-Seidel iteration is given by 


[1 0 0 0-2 2 1 
1 1 Ofu =| 0 0 =1 Jum" 4| 3| 
L2 2 1j, 0 0 0 5 
which gives 
1 0 O};J0-2 2 1 0 0j|1 
w --1 1 ojo 0O-1[u7" 4 —1 1| 013 
0-2 1/||0 0 0 0-2 I| 5 
0-2 2 1 
=|0 2-3 r+) 2 
0 0 2 -1 
The eigenvalues A of the Gauss-Seidel iteration matrix satisfy 
-à -2 2 


0 2-4 -3 |20, 
0 0 2-2 


that is 4Q — 2)? = 0, and the eigenvalues are 0, 2 and 2. (In fact, this is obvious 
since the eigenvalues of a triangular matrix are just the elements on the main 
diagonal.) Since the eigenvalue with largest modulus exceeds unity the method 
diverges. 


(b) The iterative scheme is 


0-2 2 1 
x" —|—-] 0 1 [x 4 [3]. 
-2-2 0 5 
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Taking x" = (a, B, y), we obtain 

x" = (~28 + y+ 1,-a—y 43, 2a — ff +5), 

x? =(—2a — Aff + 2y + 5,20 + 4B ~ 2y — 3,20 + 48 — 2y — 3) 
and — x? = (1, 1,1). 


The third iterate is independent of («, B, y), and x") is the correct solution of the 
linear equations. 


This result is to be expected, since p(J) = 0 so that the correction vector is 
eventually zero, i.e. the scheme terminates after a finite number of steps. 


Solution to SAQ 12 


The form of the matrix for the natural ordering is 


1 2 3 4 5 6 7 8 9 


x x x ] 
X X x x 
X x x 
x X x x 
x X X x xf, 
x x x x 
x, x x 
x X x xX 
L x x xd 


and the partitioning shows that it is not the shape of either of the required forms, 
since the diagonal blocks are not of diagonal form. 


Solution to SAQ 13 


(a) The statements are equivalent because the SOR method merely changes the 


(b) 


Gauss-Seidel adjustments by a constant factor o. Alternatively, the SOR 
iteration is given by 


(I — oL)u"*? = (1 — c)I + oU]u" + cb 
and the Gauss-Seidel by 
(I — L"*? = Uu" +b. 


In each case the points with superscript (n + 1) are related to the matrix on the 
left of the equation, and the nonzero elements in J — wLare in the same positions 
as those in I — L. 


For the natural ordering we have the following table (corresponding to the 
table in the text relevant to the diagonal ordering): 


points with superscript (n 1) points with superscript (n) 

r-l none 2,4 

2 1 5,3 

3 2 

4 1 42,9 

5 42 8,6 

6 5,3 9 

7 8 

8 75 9 

9 8,6 none 


This is the same as the table in the text. Therefore, the natural ordering 1s 
consistent with the diagonal ordering of Equation (6). 


(c) The matrix for the ordering 1, 2, 3; 6, 5, 4; 7, 8, 9 is given by 
1 2 3 6 5 4 7 8 9 


z 7 
x x x 
x x x x 
x x|x 
x|x x x 
x X X X x 
x x X|x 
x| x x 
x x x x 
L x x xd 


If we look at the computation of uf"* "'(.e., the fourth row) we see that it involves 
ust), uw and uf. This is not consistent with r = 6 in the table in the text for 
the diagonal scheme, which involves u{"*!” In the checker board scheme we 
find that ué"*” is computed from a knowledge of uf"* U, u£* ! and uf/* , so 
that there is no consistency here either. 


Solution to SAQ 14 


(a) For the problem of SAQ 10 we found [p(J)]? = 4. The optimum o for the 
SOR method is given by 


2 


Om TET 


p = 1.039. 
The spectral radius of the SOR matrix is then Wop — 1 = 0.039, considerably 
smaller than the 0.125 of the Gauss-Seidel method. 


The matrix of the equations has Property A since it has the form of Equation (8) 
of Section 11.4.3; all the block matrices are 1 x 1 matrices. It is also consistent 
since no interchanges are required to transform the matrix to this form. 


(b) The following is a sample computer output from the library program $SOR321 
in which both the Gauss-Seidel method (over-relaxation factor = 1.0) and 
SOR method (over-relaxation factor = 1.039) have been used to solve the given 
system starting with the initial approximation (u,, u3, u3) = (0, 0, 0). 


THIS PROGRAM USES THE METHOD OF SUCCESSIVE OVER RE 
" ELAX D 
TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO SIX eatin zu 


TYPE NUM3ER OF EQUATIONS(<=10)?3 


TYPE ELEMENTS OF COEFFICIENT MATRIX,ROW BY ROW 
ROW } 71,9.25,0 

ROW 2 ?20.25,15,0.25 

ROW 3 ?70,0.25,1 


TYPE ReH.Se CONSTANTS 
7105939305 


TYPE INITIAL APPROXIMATION 
70,925,0 


TYPE OVER-RELAXATION FACTOR?1.53 
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SOLUTION 


NUMBER OF ITERATIONS NEEDED IS 10 


REPEAT WITH ANOTHER OVER-RELAXATION FACTOR/INITIAL APPROX. 
YES OR NO? YES 


TYPE INITIAL APPROXIMATION 
70,0,0 


TYPE OVER-RELAXATION FACTOR?1.033 


SOLUTION 


NUMBER OF ITERATIONS NEEDED IS 7 


REPEAT WITH ANOTHER OVER-RELAXATION FACTOR/INITIAL APPROX. 
YES OR NO?NO 
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11.7 APPENDIX (Optional) 
The Eigenvalues of a Block Tridiagonal Matrix 


We show how to derive the eigenvalues of the block tridiagonal matrix A of order 
(p — 1)(q — 1) given by 


B I 
I BI 
I BI 


where B is a square matrix of order q — 1 and I is the identity matrix of order q-1. 
If A is an eigenvalue of A and x = (x1,,X12,... Xii X21 Xp- 1-1) iS a COT- 
responding eigenvector then 


[B-A I [x | 
p Bear I x2 
I B-AI X3 
=0, 
IB-a I Xy-2 
I B-Ó| xa 
where x, = Gui Xia. Xy, 1). This yields the vector recurrence relation 
Xpa + (BR ADK, +X, = 0 
together with the boundary conditions 
Xo 7 X, = 0. 
We now write v = 4 — jt where p is an eigenvalue of B with a corresponding eigen- 


vector b and suppose that 
X, = ab. 
Hence 
(0,41 — vo + t, ,)b = 0 
and, since b # 0 by the definition of an eigenvector, we have 
kti = YO + k-i = 0, 
% =a, = 0, 


or equivalently, 


-y 1 a 

] -y 1 ay 
= 0. 

1 ~y 1 Opa 

1 =y pl 
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That is, v is an eigenvalue of the Square matrix of order p — 1 


0 1 
101 
101 


L 1 0 
and is given by 


vi = 2cos = i-21L2,...,p—-1 
p 
with an associated eigenvector given by 
. in, 2i i = Ur 
sin—, sin Z, vod sn E Di è 
p p p 


as shown in the Appendix to Unit 8. 


Hence we obtain the result that the eigenvalues of the block tridiagonal matrix A 


are given by 
j= Hj vi i=1,2,...,.p—1, 
where p; are the eigenvalues of B and 


u= 2eos= i= 1,2,...,p- 1. 


The corresponding eigenvector is given by 


2in . (p — Din 


J=1,2,...,.9-1, 


- . dm ; 
x^? = | sin — bj;sin — bj; ;sin 
P p 
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Set Books 


G. D. Smith, Numerical Solution of Partial Differential Equations (Oxford, 1971). 
H. F. Weinberger, A First Course in Partial Differential Equations (Xerox, 1965). 


It is essential to have these books; the course is based on them and will not make sense 
without them. They are referred to in the text as S and W respectively. 


Unit 13 is based on W: Chapter VII, Sections 36 to 38. 


Conventions 


Before working through this text make sure you have read A Guide to the Course: 
Partial Differential Equations of Applied Mathematics. References to Open University 
courses in mathematics take the form: 
Unit M100 13, Integration II for the Mathematics Foundation Course, 
Unit M201 23, The Wave Equation for the Linear Mathematics Course. 


PDE 13.0 


13.0 INTRODUCTION 


This unit is concerned with properties of the eigenfunctions and eigenvalues of the 
system 


(py — qu + Apu - 0.— in (a, f), (1) 
u(x) = u(B) = 0, 


and of similar systems having different, homogeneous, boundary conditions. Such 
Systems are known as Sturm-Liouville systems and arise, for example, when the 
separation of variables technique is applied to solve boundary value problems for 
partial differential equations; so far in this course we have only looked at the simplest 
types, where p, q and p are constant functions. 


The differential equation (1) may have non-constant functions as coefficients and 
we may not be able to solve it explicitly using analytical methods. However, the art of 
dealing with difficult problems is to discover useful things about the solution without 
necessarily being able to express it as a formula. The present unit is a collection of 
results of this sort. Illustrations are chosen for which solutions can be written down, 


but the results apply equally to systems about which we would otherwise be wholly 
in the dark. 


The results can be summed up very briefly: the system (1) has most of the properties 
of the simplest of all such systems, namely 


u” + du=0 in (0, 1), (2) 
u(0) = u(1) = 0. 


This system has eigenvalues A4, = k?x?, and eigenfunctions u(x) = A, sin krx. The 
eigenvalues are discrete, positive and tend to infinity, and the eigenfunctions oscillate, 
u, having k — 1 zeros in (0, 1). The eigenfunctions can be used in infinite series (Fourier 
Sine series) to represent more or less arbitrary functions on [0, 1], and for sufficiently 
well-behaved functions the series converge uniformly on (0, 1]. These properties are 
all shared with the solutions of the general system (1). It can even be shown that, for 
k large, the eigenvalues of (1) are nearly proportional to those of (2). Practical instances 
of this theory appear in Unit 14, Bessel Functions. 


In order to prove these results, the system (1) will be reformulated in two different 
ways; as an integral equation in which the Green's function is used, and as a minimiza- 
tion problem. The second idea will be completely new to you and you should take a 
lot of trouble to make sure you understand it. 


Some of the SAQs have been marked as being difficult. Do not spend a lot of time 
trying to do them yourself, and if in difficulty go straight to the answer. 
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13.4 EIGENVALUES AND EIGENFUNCTIONS 


Notes 


(i) W: page 160, lines 1 to 4 
See, for example, W: Section 14. 


(ii) W: page 162, Equation (36.7) 
The system (36.1), (36.2) may be written as 
(pw) — qu = —Apu in (0, 1), (1) 
u(0) = u(1) = 0. 


Suppose the system has a solution, that is, 2 is an eigenvalue. Since all the eigen- 
values are positive (W: page 161, line — 11), 4 = 0 is not an eigenvalue and the 
system 


(py —qu 20 — in(01) (2) 
u(0) = u(1) = 0 


has only the trivial solution u = 0. Therefore, a Green's function G(x, é) exists 
for this system (Unit 9, Green’s Functions I, SAQ 9). Now, for the system (36.6) the 
equation in W: page 162, line 2 holds for every function F. It therefore holds when 


F(x) = Ap(x)u(x), 


À being an eigenvalue and u a corresponding eigenfunction of (1). Thus we have 
(36.7). 


Equation (36.7) is an integral equation, which means that the unknown function 
occurs under an integral sign. Neither u nor A is known in this equation. Note 
that the boundary conditions are automatically incorporated in (36.7), since the 
integral tends to zero as x ~0 or x ~ 1. 


In the subsequent work, you should notice that we do not have to know what G is, 
but only that there exists a Green's function belonging to the system (2) and 
that it has some simple properties. 


(ii) W: page 162, line — 10 
We met the result that 


E 
> a, converges = lim a, = 0 
i nao 


in Unit M201 20, Convergence and Bases. Here we put a, = 1/22, and deduce 
that 42 grows unboundedly. Since 2, > 0, we have 4, —> cc. 


It is assumed that the eigenvalues of the given system form a discrete set of 
numbers. This will follow from the minimum principle described in Section 13.1.2. 


SAQ I 

Find the eigenvalues and eigenfunctions of the system 
uw +iu=0  in(wf) 
u(x) = u(B) = 0. 

(Solution on p. 28.) 

SAQ 2 

Find the eigenvalues and eigenfunctions of the system 
u" +2u=0 in (a, p), 
u'(a) = u'(B) = 0. 
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In W: page 161 it was proved that the system (36.1), (36.2) has no eigenvalue 4 = 0, 
but 4 = 0 is an eigenvalue of the system above. Where does the theory need changing 
for this system? 


(Solution on p. 28.) 


SAQ 3 $ 
The system 

(pé) — qu + Apu-0 — in(a, f), 

u(a) = 0, u(f) = 0, 


has an eigenvalue 4, and a corresponding eigenfunction u,. Show that every eigen- 
function of 4, is a multiple of u,, that is to say, the dimension of the eigenspace is one. 


HINT: Suppose there are two linearly independent eigenfunctions. Then every solution 
of the equation with 4 = 4, must be a linear combination of these. This leads to a 
contradiction. B 


(Solution on p. 29.) 


SAQ 4 


For the following system, prove that the eigenfunctions corresponding to different 
eigenvalues are orthogonal with respect to the weight function p. 


(pu) — qu + Apu = 0 in (a, f), 
Au'(x) + Bu(x) = Cu'(B) + Du(f) = 0, 
where A and C are nonzero. 


(Solution on p. 30.) 


SAQ 5 
Find the eigenvalues and eigenfunctions of the system 
du . 

4p ^9 —n<O0<n, 


with periodic boundary conditions, 
u(—n) = u(n), u'(—2) = u(x). 


Show that the solution space corresponding to each nonzero eigenvalue has dimension 
2. 


(Solution on p. 30.) 


SAQ 6 
Given the system 
uw. 4 àu-0 in (0, 27), 
u(0) = 0, 
u'(2x) — u(2x) = 0, 
show that the eigenvalues are the solutions 4 of the equation 
At = tan(2n43), 
and that the least eigenvalue is negative. 


(Solution on p. 31.) 
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13.1.2. The Minimum Principle 


In Section 13.1.1 we succeeded in obtaining certain properties of the cigenvalues by 
replacing the differential equation system (36.1), (36.2) in W by an integral equation. 
In this section we reformulate the problem in another way which may be completely 
new to you. It is a method of wide application, and represents a special case in the 
subject called the calculus of variations. The simplest problem in this subject is to 
find a function which minimizes the numierical value of a given definite integral in 
which the function appears. For example, the function minimizing 


x/2 
S WO- 92601 ax 
o 


( 


X) = sin x. 


peu 


Notes 


(i) 


(ai) 


(iii) 


W: page 163, line 6 

The passage leading up to (36.8) is intended merely to motivate an approach 
which begins with (36.8), since suddenly to begin examining the form (36.8) with- 
out any suggestion that it might be profitable to do so would be unreasonable. 
There is therefore no need to learn this argument. 

W: page 163, lines 8 to 10 

This is a usc of Parseval’s equation in its general form, W: page 75, Equation (16.6), 
with 


loe - adj 


in place of f *. 


Since 4, > a, for n = 1,2,3,..., we now have 


wo 1 00 i 1 3 

4 4» 2 2 Z 1 
Daaf manfa pu = 2, f po’, 
n=1 o n=1 0 0 


by the equation on line 7. 
W: page 163, line —12 


A function is said to be piecewise continuously differentiable or piecewise smooth 
on a finite interval [«, £] if the interval can be divided into a finite number of 
open subintervals (x, x), (x1, x2), .... (x,.,, f), such that the function and its 
first derivative are continuous on each subinterval and the right- and left-hand 
limits of the function and its first derivative exist at the end points of each sub- 
interval. This class of functions embraces most functions directly encountered 
in physical applications. The following figure shows a continuous, piecewise 
continuously differentiable function. 


(iv) 


(v) 
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W: page 163, line — 10 
For example, put P=1,q¢=0,p=1.The system becomes 
u” + hu=0 in (0, 1), u(0) = u(1) = 0, 
with eigenvalues A = n?n? 


is suggested (not proved) i 
piecewise continuously di 


Pe 
= >n’, 
0 


For example, if $(x) = x(1 — x) we obtain for the left-hand side 


» the first eigenvalue 2, being equal to 22, Then what 
n the preceding paragraph is that for any continuous, 
fferentiable function 9 satisfying $(0) = (1) = 0, 


1 
Í (1 — 4x + 4x?) dx 
0 


1 
I (x? — 2x3 + x4)dx 39 
0 


which is greater than z? ~ 9.87 by a small amount. You may try other functions, 
but will always fail to produce a value less than 72, although the function 
$(x) = sin x gives equality. Therefore the result seems likely to be true. 


So, instead of trying to refine the given proof so that the gap in line 6 is filled, we 
restart from the standpoint of Equation (36.8), by proving directly that if there 
exists an admissible function $ which minimizesthe value of this expression, then 
the minimum value attained is the lowest eigenvalue, and the function $isa 
corresponding eigenfunction. 


The proof that there is a minimum, given by some admissible function, is beyond 
the scope of this course. Generally speaking, to prove the existence of something 
without being able to display it is a difficult matter. It is easy to show similar 
expressions that have no minimum: for instance, 


1 
Í $? 
o 
i 
fo 
0 
has no minimum. We shall simply assume throughout this unit that there does 
exist an admissible minimizing function for the Rayleigh quotient. 


W: page 163, line —7 


$ is admissible if à is continuous on [0, 1], $' is piecewise continuous on (0, 1], 
and $(0) = $(1) = 0. 


W: page 163, lines —5 to —3 

V gives the minimum value, y, to the Rayleigh quotient, and therefore for any 
other admissible function, such as y + to, the “2” in line —5 holds. (Note that $, 
though still an arbitrary admissible function, has a different significance from 
its earlier use.) In particular, for any given admissible ¢, the “2” holds for every 
value of z, and becomes an equality when t = 0. Therefore, if ¢ is a given function, 
the left-hand side is a smooth function of t with a minimum at t = 0. 


(vii) W: page 164, lines 2 to 4 


The integration by parts gives 
H 1 
- [| ettpY — ay + nol + [ove] =0. 
o 


The bracket vanishes because (0) = (1) = 0; this is the stage in the argument 
where the boundary conditions appear. Note that we are now assuming that a 
minimizing function y not only exists, but also is twice continuously differentiable 


throughout (0, 1). 
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(viii) W: page 164, last paragraph . 
Let $ = w be that function which minimizes the Rayleigh quotient amongst all 
continuous, piecewise continuously differentiable functions $ satisfying 


$(0) = $(1) = 0 


and 


1 
MED (a) 
o 
Now the proof proceeds exactly as before up to W: page 164, line 12, i.e. 
(PWY — qi + upy = 0, 
y0) = Wl) = 0. 


is therefore an eigenfunction which is, by (1), orthogonal to u,. It must there- 
fore be an eigenfunction with an eigenvalue different from 4, (by SAQ 3) and 
the smallest available eigenvalue is 2,, with eigenfunction uz. 


The proof for subsequent eigenvalues is treated in the same way. 

(ix) W: page 165, lines 4 to 7 
This shows that, for any k, the family of admissible functions which are orthogonal 
to the first k — 1 eigenfunctions is nonempty—a polynomial is constructed. That 


is to say, no matter how large k may be, we have admissible functions for the 
next step. 


SAQ 7 
Verify by substitution that 


where A, is the lowest eigenvalue and u, any corresponding eigenfunction of the system 
u" + du =0 in (0, 1), 
u(0) = u(1) =0. 

(Solution on p. 31.) l 


SAQ 8 

Show that the smallest eigenvalue of the system 
u” + àu = 0 in (0, 1), 
u(0) = w'(1) — u(1) = 0, 

is arrived at by minimizing 


f $? - 0 


o 
1 
Le 
0 
over all continuous piecewise continuously differentiable functions $ satisfying 


(0) = ¢'(1) — (1) = 0. You may assume that a twice continuously differentiable 
minimizing function exists. 


Ri) = 


(Solution on p. 32.) 
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SAQ 9 


Let A, be the kth eigenvalue and uy, 


the kth ei i i iodi 
boundary conditions, eigenfunction of the system with periodic 


u" + Au =0 in (—z, x) 
u(~z) = u(x), u'(— n) = u'(n), 
which you solved in SAQ 5. (Double eigenvalues are to be counted twice.) 


Proceed as in W: pages 163 to 165 to show that 
= $9" 
A, = min — 


p? 


the minimum being taken over functions satisfying 
$(~n) —$(x), $'(—n) = ó'(x), 


and subject to the orthogonality conditions 


fon =.. = J gun =0. 


(Solution on p. 32.) 


SAQ 10 (Very difficult) 


Prove that the minimum of the quotient 


0 


over all smooth functions satisfying only the condition $(0) — 0, is equal to the 
lowest eigenvalue of the system 


(pu) — qu + Apu = 0 in (0, 1), 
u(0) = v'(1) = 0, 


and the minimizing function is a corresponding eigenfunction. By considering a 
particular choice of p, q and p deduce that 


1 1 
Í o? 2 in? f p, 
0 0 


for every smooth $ for which $(0) — 0. 


(Solution on p. 33.) 


SAQ 11 
A certain wave propagation problem inside a sphere of radius a gives rise to the 
following system: 


du 2du 
dé gs bw O0<r<a, 


u(a) = 0, u(0) is finite, 
where r is the distance from the centre. 


(a) Find all the solutions. 
HINT: Determine the differential equation satisfied by ru. 


. M 


(b) Suggest a minimization problem equivalent to this onc. (Remember to use the 
self-adjoint form.) 


(Solution on p. 34.) 
SAQ 12 (Difficult) 


Let p, q and p have the same meanings as in W, and let v, w be any two piecewise 
continuously differentiable functions on [0, 1]. Now define E and H by 


1 
E(v, w) = f (po'w + quw), 
o 


and 


1 


H(v,w) = Í pow. 
0 


(a) Write the Rayleigh quotient in terms of E and H. 
(b) Prove that H defines an inner product.* 
(c) Prove that E defines an inner product when q(x) > 0 for all x e (0, 1). 


(d) Carry out the proof that the minimum of the Rayleigh quotient is equal to the 
lowest eigenvalue (W; page 164, Equation (36.10)) using this notation. 


(e) Show that any function solving Problem B also solves Problem A: 
Problem A 
Minimize 
1 
li (pp? + q$?) 
2 
f po? 
0 
over all smooth functions for which $(0) = $(1) = 0. 
Problem B 
Minimize 
1 
f Gor + 40% 
0 
over all smooth functions for which $(0) = $(1) = 0, and 


foe -1 


0 
(Solution on p. 34.) 


* The definition of an inner product is given in the Appendix to Unit 3, Elliptic and Parabolic Equatior 
aS. 
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13.2 THE CONVERGENCE OF GENERAL FOURIER SERIES 


13.2.4 The General Theory 


The detailed analysis involved in this section is very difficult, and we therefore recom- 
mend that you simply read through it, noticing the principal results. Some of the details 
are not spelled out in W, and in Section 13.2.2 we prove these results. If you are pressed 
for time you should omit this section altogether and continue with Section 13.2.3. The 
general nature of the results is as follows. Let {u,} be a sequence of eigenfunctions of 
the system (36.1), (36.2) in W: page 160. It has been proved in Section 13.1 that they 
form an infinite set of orthogonal functions. We have seen in Unit 6, Fourier Series 
that it is possible to expand an arbitrary function f on the interval (0, 1) in terms of 
thesc eigenfunctions. We should like to be able to say further that if 


$ 
J Pit) ax 


o 


1 
Í P(x)ug(x) dx 


o 


then 
fo =E gu) O<x<1, 
1 


and that the convergence is uniform (this being essentially the result for trigonometric 
Fourier series provided f is suitably restricted). 


In the next reading passage it is shown that the eigenfunctions {u,} are as good as 
the trigonometric functions for this purpose. In fact the convergence properties of 
trigonometric Fourier series follow as a particular case of the more general series 
above. 


Notes 
(i) W: page 165, line 16 
The functions u; and u; are orthogonal when i gy 


(i) W: page 165, Equation (36.12) 
Note that, for example, 


k—i 2 k-1 k=1 
5 ct =% an Y dd 
1 n=1 m=1 

k—1 k-1 


Y M eeu, 


n=1 m=1 


(ii) W: page 166, line 2 
Since the terms in the series are positive the right-hand side is bounded above by 


d qr. 
zl (nf? + af) 


which approaches zero as A, ~ oo. 


(iv) W: page 166, lines 6 to 8 
The proof that any square integrable function can be approximated in the mean 
arbitrarily closely by a continuously differentiable function is beyond the scope 
of this course. (However, if you have ample time at your disposal, you may like 
to sketch out a proof that a piecewise smooth function can be approximated in the 
mean arbitrarily closely by a continuously differentiable function.) The proof of 
completeness is given as Theorem A in Section 13.2.2. 
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(v) W: page 166, lines 9 and 10 
This statement means that if 


H 
MA 
0 


Cn = M n= 1,2,3, 


2 
MN 
o 


t k 2 
lim e[s- » cat] =0. 
kom Jo r= 


(vi) W: page 167, lines 1 to 3 
From W: page 166, line —4, with f(£) = G(x, ë), we have 


then 


oo 
G(x, č) ~ È eu 
1 
where the Fourier coefficients {c,} depend on x. They are given by 


1 
| (6t cer dt 
2:0 


€, 


f POE dë 


according to the general formula at the bottom of W: page 162. 
But by Equation (36.7), with A = 4, and u = u,, 


1 1 
[, POO uE dě = 1-9. 
o n 
and therefore 
1 u,(x) 
1 X9 0 
fp tule? dc 
0 


C, = 


as in line 2; d, is treated similarly. 


(vii) W: page 167, Equation (36.14) 
Here we have expressed the Green's function for the system (36.6) in terms of the 
eigenvalues and eigenfunctions of an associated eigenvalue problem (36.1), (36.2), 
in which p is any positive continuous function of our choice with domain (0, 1]. 
This formula is called the bilinear formula for the Green's function. Note that the 
symmetry of G is clearly displayed in the series. 


(viii) W: page 167, line —9 
Schwarz's Inequality in a finite dimensional Euclidean space gives 
k 2 k k 
PORPP 
i21 iz1 i=1 


for any real numbers a,,..., âp, bis- .. , dy. 


General Comment 
From SAQ 9 we see that the trigonometric functions, 
cos nÜ, sin nÜ meas (1) 


are the eigenfunctions of a system of the type considered, but with periodic boundary 
conditions, on the interval [— x, n]. They also satisfy a minimization principle similar 
to (36.11), with the appropriate modification to the boundary conditions. By making 
the small changes in the theory arising from the changes in the interval and boundary 
conditions we can show that the set of functions (1) is complete on [ — z, z), and also 


deduce the standard theorems on convergence which we quoted without proof in 
Unit 6. 
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SAQ 13 
The Green's function for the problem 
u"(x) = — f(x) Ozxzcl, 
is given by 
1-4 
6x8 = f T 
g&l- x) xz 


(Example, W: page 122 
series on [0, 1]. 


[: 
[s 
). Use the bilinear formula (36.14) to express G as a Fourier 


HINT: Make a suitable choice for p. 


(Solution on p. 35.) 


SAQ 14 


Prove that the set of functions sin nx (n = 1,2,...) is complete on 0 € x «€ 1. (This 
is essentially the set used for expansions in sine series; see W: Section 20. The con- 
vergence properties used there follow from the results of this section.) 


(Solution on p. 36.) 


13.2.2 Three Theorems (Optional) 


We now fill in the details of three main results which were stated in W with only brief 
proofs. You should not spend time on this section at the expense of tlie remainder of the 
unit. 


The three results are: 
A the eigenfunctions {u,} of the problem (36.1), (36.2) are complete (W: page 166); 


B the Fourier series of a continuous function f with f(0) = f(1) = 0 and 


1 
Í, (pf? + qf?) 


finite, in terms of the complete set of eigenfunctions {u,}, converges uniformly to 
J on [0, 1] (W: page 167); 


C the solution of the nonhomogeneous problem is given by the formula in W: 
page 168, line 5. 


We first quote the Triangle Inequality* for the space of functions square integrable 
on (0, 1) with weight function p; this tells us that for any two such functions g and h 


we have 
1 + 1 H 1 + 
2 < 2 hey. 
li e + n | sili pe «lf pr | 


This result will be used in the proof of Theorem A. The essential point to be proved 
in this theorem is that a function f can be approximated in the mean by its Fourier 
series in terms of the eigenfunctions (u,). It has already been shown in W: page 166, 
line 4 that if f is continuously differentiable then this result holds. We wish to extend 
the result to more general functions which arise and can themselves be approximated 
in the mean by continuously differentiable functions. 


*The Triangle Inequality may be deduced from Schwarz's Inequality in any inner product space: for 
example, D. L. Kreider er al.. An Introduction to Linear Analysis (Addison-Wesley, 1966), p. 264. 
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THEOREM A 
1 
For all f such that | pf? is finite, 
0 


1 k 2 
lim f olr- Feun) =0, 


o 
where {c,} is the sequence of Fourier coefficients of f with respect to {u,}. 
Proof 


Let f be any function of the class considered. Then given any e > 0, let f be a con- 
tinuously differentiable function such that 


1 
J| - Pes and fo) = fu) = 0. a) 
0 


(As-stated in note (iv) of Section 13.2.1 the proof that such an f exists is omitted from 
the course.) Now let (d,) be the sequence of Fourier coefficients for f, and put 


k k 
$6) = Xon) — o0) = X dux). 
1 1 


for the kth partial sums. For any given k the integral 
1 k 2 
[ole xn] 
0 1 


is minimized when the sequence (b,) is the sequence of Fourier coefficients (W: page 
71). Therefore, for every k, 


1 1 
o« [ t «ou ar. Q) 


The last integral can be written as 
l 
Í, af-f+i-oay, 


so that 
p1 


0s Js pf - s 


1 
S| psf-ftf- a) 


vo 


Y own Feat ae jm 
< li p - fy ] + [i pu - aa] } by the Triangle Inequality. 


Therefore, taking the limit as k ~ oo, 
i 


1 
o« jm [otf sr f -a 
0 o 


by (1), where we have used the fact that 
1 
lin | p(f— a)? =0, 
ko 0 
by W: page 166, lines 3 and 4. Since e is arbitrary, we must have 
1 


lim | p(f — s? =0, 


kao Jo 


as required. 
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THEOREM B 


If f is continuous on (0, 1), (0) =f(l)=0 


H 
f S? + qf?) 


is finite, then the Fourier series 


E] 
Yu, 
1 


converges uniformly to f. 


and 


Proof 


Define the partial sums 


sx) = You) 
and 
k 1 
o= ES Í pu. (3) 
Then Equation (36.15) becomes, after taking the square root of both sides, for all 
M > 0 and for all N > M, i 
Sy(x) — suG9I < (G(x, x)| loy — oy. (4) 
Now G(x, y) is bounded; in particular, there exists a number A such that 
G(x, x) < A O<x€<l. (5) 


Also the series (3) converges as k ~ co (W: page 166), and so given any number c > 0, 
there exists No independent of x such that 

e? 

A 


Equations (5) and (6) together with (4) give 


sa(X) — Sy(x)| < S) =£ 


On — Oyl < whenever N > M > No. (6) 


whenever 
Osx<1 and N> M » Ne. 


Therefore {s,(x)} is a Cauchy sequence for every x and so converges on (0, 1] to S(x), 
say. The convergence is uniform, by Cauchy's Test (Unit 6) since No is independent 
of x on [0, 1]. Therefore S is continuous, by the properties of uniform convergence, 
since each s, is continuous. There remains the question of whether the function S 
which is the sum of the series is in fact equal to f. We know (Theorem A) that the 
series {s,} converges in the mean to f, and we shall prove the following: 


given that f is continuous, 
1 
jim p(f - s =0 (7) 
TN 


and 


jim s(x) = S(x) uniformly on [0, 1], (8) 


then 
S(x) =f(x) — xe(0 1) 


PDE 13.2.2 


To prove this we consider the integral 


1 
f s-s- f pif = s) +r- SP 
o 


where S is continuous and therefore integrable. Hence, by the Triangle Inequality, 


o< fi olf - S? «I ZEE + [fons] 


By conditions (7) and (8) the right-hand side tends to zero as k ~ co. Therefore 


1 
| os- s= 
0 
and so 
J= =S(x) O<x<l. 


For suppose that f(x) + S(x) at some point x. Then by the continuity of f and S, 
p(f — SP > Oin ne intervat containing x, and p(f — S}? > 0 elsewhere, so that 


1 
f o-s >o 
o 
contradicting the equality above. Therefore f(x) = S(x) at every point on [0, 1]. 


COROLLARY 


oo 


» c,uj(x) is absolutely convergent, and moreover » le, I (x)] is uniformly convergent. 


Proof 


The absolute convergence of the series 


Xon) 
1 


can be proved on the same lines. Simply notice that (36.15) is still valid if 


N 2 
| Y. ded bcs 


M*1 


is substituted for the left-hand side. It follows that 
D lend lnl 
1 


is uniformly convergent. (This result is not true unless f(0) = f(1) = 0.) 


These results are needed to prove the following result (W: page 168, line 5). 
THEOREM C 

1 ^" 
Let Í |F| be finite Then the problem 

0 


(pwY — qw = —F, 
with 
w(0) = w(1)=0 
has the solution 
w(x) = DET 
1 A, 
where u, and 2, are the eigenfunctions and eigenvalues of the problem 
(pu) — qu + Apu = 0, 
u(0) = u(1) = 0, 
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for some continuous positive function p on (0, 1], 


F(x) 2 : 
^ È euo. 


and the coefficients {c,} are defined by 


The series for w converges absolutely and uniformly. (Note that many expansions 
of the solution are possible depending on the choice of p) 


Proof 


Let G be the Green's function for the System; then the solution of our problem is 
given by 


1 
wo) = S Ges OF ae, 
o 


Using the bilinear formula for G we have 


ws) = [ rO ue] d 6) 
0 1 


0 


for each x e [0, 1]. By Theorem B, the Biven expansion converges to G(x, č) uniformly 
on 0 € č < 1, since the function ë — G(x, č) is continuous, piecewise smooth (so 
that it satisfies the integral condition in Theorem B) and zero at the end points 
(W: page 122). It is known that in this case we can interchange the integral and the 
summation in (9), giving 


t| pu 
0 
e ud) [ EO o 
ee otal ud) dé 
- 7 zx pu? [ p) 
0 
- LL u(x), 


according to the definition of c,. It is clear that this expression is the Fourier series 
for w. Now, w is continuous and vanishes at the points 0, 1 and 


Il 


1 
-Í w[(pw) — qw] integrating by parts 
o 


1 
1 wF, 
0 


which is finite since w is continuous on the closed interval (0, 1] and 


1 
Í IF| 
o 


is finite. Hence, by Theorem B and its corollary, the series for w converges absolutely 
and uniformly. 


1 
f (pw? + qw?) 
0 
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13.2.3 Vibration of a Variable String 


Tt is rather hard to find equations of the form (36.1) which have non-constant coef- 
ficients but which can be solved in terms of elementary functions to display the eigen- 
functions and eigenvalues we have been talking about. We shall now meet an example 
which, though artificial, is one that can be worked out in full. In Unit 14, Bessel 
Functions another example is given which requires us to explore completely new 
functions. 


Notes 


(i) W: page 169, Equation (37.2) 
Here p = 1, q = 0, p(x)-(1-?. 


(ii) W: page 169, lines — 12 and — 11 
We test whether there is a value of a satisfying the equation by putting 


X(x) = (1 + x) 
into (37.2). Then a must satisfy 
a(a — 1)(1 + xj? + 2( + x^? = 0, 
for all x. This is equivalent to 
a(a — 1) 3-4 — 0. 
(iii) W: page 170, line 5 
We have 
(1 + xj TVIRD L4 4 xg xvii, 
where $4 -—1 is real. Now, a^, when a and b are real, is defined in terms of 
the complex exponential function as e”!°**, Here a = 1 + x, b = i4 = 1. 
(iv) W: page 170, line 8 
Formally, we write the complex solution given in lines 6 and 7 in the form 
X = X, + iX,, where 
X(x) = (1 + x cos(\/(4 — Hlog(1 + x), 
X(x) = (1 + x} sin(/(4 — 3)log(1 + x), 
X, and X, being real functions. Then, since X is a solution, 


a 
ee, 
*tüxxy 


À A 
=|X%~)+—" x il X5(x) + —— f 
| T(x) + a+ i| + | 5x) + +x E 
Each of the brackets on the right is real, and therefore the complex number on 
the right is zero if and only if the brackets are separately zero. Therefore X, 
and X, are, separately, solutions of the original differential equation. Strictly 
speaking, the correct thing to do at this stage would be to verify that the real 


functions X, and X; are two independent solutions of the differential equation 
(37.2). 


0 = X"(x) (x) 
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13.3 FURTHER PROPERTIES OF EIGENVALUES 
AND EIGENFUNCTIONS 


13.3.0 Introduction 


In this section we shall mect a number of important theorems regarding the eigen- 
values of boundary value problems and the zeros of their eigenfunctions. The two 
results concerning the cigenfunctions are the Separation Theorem and the Oscillation 
Theorem, and the result about the eigenvalues is the Monotonicity Theorem. 


We begin the discussion by proving a comparison result which will turn out eventually 
to be a particular case of the Monotonicity Theorem. We then discover that the first 
eigenfunction of a boundary value problem has no zeros between the end points, 
i.e. it does not change sign on the interval. Finally we plough through the main 
theorems. 


13.3.1 A Comparison Theorem for Sturm—Liouville Problems 


We have not asked you to read the first two paragraphs of W: Section 38 because the 
argument is difficult to follow. We shall however obtain Weinberger’s result in 
Section 13.3.2. 


Notes 


(i) W: page 172, line —4 
Note the shape of the trial function, ¢. 


A 
dx) 


(ii) W: page 173, lines 3 to 5 
The ratio referred to is 


B 
(pb? + qd?) 


i pd? 


The reasoning is not easy to follow, but is credible. The strict inequality 
Ay S Ay 


is important. 
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General Comment 


The two arguments which we have described as difficult to follow (W: page 172, line 4 
and W: page 173, line 3) are valid, but they rely on a slight extension of the result in 
W: page 163. Effectively, it can be shown that if à is continuous and $" is piecewise 
continous on [«, fj], 


(a) = 96) 2 = b(x,) = (8) = 0 


where X,,...,X,, are the points of discontinuity of d" and the Rayleigh quotient 
R[@] = A, then ¢ is an eigenfunction corresponding to /,. 


13.3.2 The Zeros of Eigenfunctions 
In W: page 172 (the first paragraph), it is shown that the first eigenfunction, u,, does 
not vanish on (x, f). To express this another way, 
the first eigenfunction of the system 
(pu) — qu + Apu =0 in (o, ff), 
u(x) = u(f) = 0, 
is of one sign throughout the interval. 
Of course, we can change this sign by multiplying the function by a constant. 


We shall see in the proof of the Oscillation Theorem that the kth eigenfunction has 
at most k — 1 zeros; thus our result is just a special case of this. 


Tt follows that 
any other eigenfunction changes sign on (x, B), 
Li 
for if k > 1, Í piu, = 0; or alternatively 


if a and f| are consecutive zeros of any solution of the differential equation 


(pu — qu + Apu = 0 


with any given 2, 4 must be the lowest eigenvalue for the interval [a, p]. 


Notes 


(i) W: page 173, SEPARATION THEOREM 


The proof given is not complete, since it assumes that u has at least two Zeros. 
Substitute the following proof. 


Proof 


The proof is by contradiction. Suppose that u does not vanish on the open 
interval (a, J); then u and v each has constant sign on (a, 8) and we may assume 
without loss of generality that u > 0, v > 0 on (a, B). Multiply Equation (38.4) 
by v and Equation (38.5) by u and subtract, then integrate from & to f (compare 
Lagrange's Identity, Unit 9). We obtain 

e B 

], POY —upey <2 — a [ pw <0, () 


with strict inequality if A < 4. 


Now integrate the left-hand side by parts. We obtain 


p 
L — woe) <0, 


z 
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and since v(@) = vB) = 0 by hypothesis, this becomes 


p(2)u(&)v(a) — p(B)u(B)v'(B) < 0. 


Since v > 0 on (&, B) and v(&) = v(B) = 0, we have v'() > 0 and v'(B) < 0. Now 
suppose firstly that u(3) and u(B) are not both zero. Since p(x) > 0 on [, B] 
(W: page 160, line 12), and u(x) > 0 on [à, J], the left-hand side of (2) is positive. 


(2) is therefore contradicted, so u must change sign, and therefore vanish, some- 
where on (à, ). 


(2) 


Secondly, suppose that u(@) = u(B) = 0 and 2 Æ 2 (in this case u and v are 
independent eigenfunctions for [&, B]. Then (2) becomes 

Pua) — p(Byu(B)v'(B) < 0, (3) 
(see Equation (1), and again there is a contradiction, since the left-hand side is 
zero. 


Finally. suppose that u(@) = u(B) = 0,and 4 = 7. Then wand vare eigenfunctions 


for the interval [&, 7], with the same eigenvalue, and are therefore multiples of 
each other (see SAQ 3). 


W: page 174, lines 1 to 9 
We expand the argument in this note. 


Firstly, let {c,,} be any set of constants, not necessarily satisfying the orthogonality 
conditions on W: page 173, lines —4 and, —3. Then, 


p p l 2 I 2 
[oen - PLE Sm] Ee] 


m-l 
B 1 d i 1 
=| PE D 007 4 +g YY pepu | 
a m=l n=l metal 
Since uu = Wy = 0 for m # n, this is equal to 
A t à $ SICURA 
27, 07" 3 5 
ll L Y (um)? +4 Y, ciu? |. 
a mz1 m=) 


and since u"(x) = 0 except on (X,,-1,X»,) where it equals u,(x), this is in turn 
equal to ` 


Xm 


1 
> Ge [pug + qui). 


m=i Xmas 


The argument continues up to line 4. Therefore, finally, 
B 
Í (pb? + qd?) 


Now, (1) holds for any set of constants c, (m - 12,.... 1) We can shere ore 
choose these | constants to satisfy the | — 1 linear equations W: page 173, 
lines —4 and —3, 


B " me 
[ ot =0. f pdms = 0: a 


that is, we can arrange for $ to be orthogonal to the first l-1 UMS 
without disturbing the truth of (1). But the absolute minimum o 


(1) 


u 2 
| (po? + qd?) 


r po? 


ns] 
[n 


over all admissible $ subject to (2) is, by the minimum principle, equal to 4j. 


Therefore 
dy S ho 


and so | < k. Now the assumed number of zeros of u, on the open interval (a, fj) 
is equal to / — 1 ( W: page 173, line — 11). Therefore the number of zeros of ug 
on (a, fl) is not more than k — 1. But we know that the number of zeros is at 
least k — 1, so that equality must hold ( W: page 173, line — 14). 


Note that, in particular, it follows that the first eigenfunction has no zeros in (a, p). 


SAQ 15 

Let up, Ugi be successive eigenfunctions of the system 
(pu) — qu + Apu = 0 in (0, 1) 
u(0) = u(1) = 0. 


Show that there is exactly one zero of u, on the open interval between every two 
zeros of t4, in (0, 1), and exactly one zero of u,,, between every two zeros of up 
(i.e. the zeros interlace). 


(Solution on p. 36.) 


SAQ 16 
Let u*, u** be any two linearly independent solutions of the equation 
(pu) — qu + Apu = 0 in (0, 1). 


Prove that on the open interval between any two successive zeros of u* there lies 
exactly one zero of u**, and on the open interval between any two successive zeros of 
u** there lies exactly one zero of u*. 


(Solution on p. 36.) 


SAQ 17 


Confirm that the nth eigenfunction, X,, of the variable string problem in W: Section 37 
has exactly n — 1 zeros in (0,1). Check that the zeros of X, and X,,, interlace 
(see SAQ 15). 


(Solution on p. 36.) 


SAQ 18 
W: page 175; Exercise J. 
(Solution on p. 36.) 


SAQ 19 (Difficult) 


(a) Prove that, ifq = 0, reducing the interval («, fj), increasing p, q or b, or decreasing p 
increases the lowest eigenvalue of the system (38.6). (You may assume that the 
minimum principle (38.7) applies.) 


(b) W: page 175, Exercise 2. 
(Solution on p. 37.) 


SAQ 20 (Difficult) 


Prove that the kth eigenfunction for the system (38.6) has at most k — 1 zeros on the 
open interval « < x < ff. 


(Solution on p. 38.) 
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SAQ 21 (Difficult) y 


Prove in full the result in W page Ines . nvalue: 
g 5. lines 16 
í ^ ( j 175.1 and 17. that the e genvalues o systern 


(Solution on p. 39.) 


SAQ 22 
W: page 176, Exercise 3. 
(The meaning of the result is that 2, behaves like k?z? for large k.) 


HINT: Compare the system with a system having constant coefficients. 
(Solution on p. 40.) 


SAQ 23 
Consider the system 
u” + Apu -0 in (o, f), 
u(a) = u(f) = 0. 
(a) Prove that 
k?n? <h< en? 
(B=) Pm © (B — e pus 
where Pmax Pmin represent the maximum and minimum values of p respectively 


on (a, fj). 


(b) Let the kth eigenfunction for the system have zeros at the points « = xo, Xi, Xs. 
oo Xy = B. Prove that for i =0,1,....k — 1, 


n? x 


——— i AM MÀ 
c "uv sek 2 
(Xii — XO Pimax (iei — X Pismin 
where p; max and P; min are the maximum and minimum values respectively of p(x) 
on the interval [x;, x;4,]. Deduce that 


> 


x 
Pimia — X) S p S plioaXie 1 — Xi 


fori20,1l...,k—- I. 


(c) Let M, be the maximum distance between successive zeros of the kth eigen- 
function; that is, the maximum value of (xj; — x). i 0 L...., k-i. 
Assuming that jim M, = 0, prove that 

E 


NE^ T 
lim => = 


a 
ko k? fi pit) ax] 


(Solution on p. 40.) 


134 SUMMARY 


This unit is concerned mainly with the properties of the eigenvalues and eigenfunctions 
of the self-adjoint system 


(pu) — qu + Apu = 0 in (a, f), d) 
u(a) = u(B) = 0, 


and with related systems having boundary conditions involving u'(a) and u'(); p,q and 
p are continuous functions on [a, b], p has a continuous first derivative on (c, ff), and p 
and p are positive and q nonnegative on [a, fj]. 


Some of the results in W are stated for the standard interval [0, 1]. Since they apply 
without change to an arbitrary interval [«, fj], we quote them for this general interval. 
The eigenvalues 4,, 25,..., are assumed to be indexed in order of size. 


THE MINIMUM PRINCIPLE 


Let $ be a continuous, piecewise smooth function on [o, fj], and let $(a) = $(fl) = 0. 
Then 


à, € 8 —,— —, Q) 


where 2, is the lowest eigenvalue of the system (1). If $ is allowed to range over all 
continuous piecewise smooth functions, then the minimum of the right-hand side 
is attained when @ = u,, where u, is an eigenfunction corresponding to 4, , and then 
equality holds in Equation (2). 


For the other eigenvalues we have the following extension: 
d 2 2 
f to? + 46%) 
4, = min 2——7,— ——. 
| pg? 
a 

the minimum being over functions $ which satisfy the conditions 
(a) œ is continuous and piecewise smooth on [c, f], 


(b) (a) = $(f) = 0, 


Li Li Li 
(c) | pou, = | pou, =- = Í pou, = 0, where u,,... , u are eigenfunctions 
a a a 
corresponding to the first k — 1 eigenvalues, 
Any minimizing function of this class which is twice continuously differentiable on 
(a, B) is an eigenfunction. 
Properties of the eigenvalues and eigenfunctions of the system (1) 


The eigenvalues are positive, infinite in number and discrete (this justifies indexing 
them in order). 


Eigenfunctions belonging to a single eigenvalue are linearly dependent. 


Eigenfunctions corresponding to different eigenvalues are orthogonal with weight 
function p, i.e. 


B 
I Puyt = 0, ksl 


MONOTONICITY THEOREM 


Reducing the interval [z, fj] to a subinterval of [a, fj], increasing p or q, or decreasing p 
increases all the eigenvalues. 


26 
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SEPARATION THEOREM 
Let u and v be any two solutions of the equations 
(pu'Y — qu + Apu = 0, 
(pv) — qv + 2p» =0 


respectively, with 4 < 4. Then if à and 


y 1 B are consecutive zeros of v, there is at least 
one zero of u in the open interval ( 


&, B), unless 7 = 2 and visa multiple of u. 


OSCILLATION THEOREM 


The kth eigenfunction u, has exactly k — 1 zeros in the open interval (o, f). The zeros 
of successive eigenfunctions interlace. 


THEOREM A 


The set of eigenfunctions {u,} is complete for the space of functions f for which [ pf? 
is finite; i.e. for such functions i 
, R k 2 
im f olr- DEM =0 


where {c,} is the sequence of Fourier coefficients of f with respect to {u,}. 


THEOREM B 
If f is a continuous function such that f(a) = f(f) = 0, and 
B 
for san 
a 
exists, its Fourier series in terms of the complete set of eigenfunctions (u,) converges 
uniformly to f on (a, fj]. 
We also considered briefly the more general problem 
(ply — qv + upo 2 0 — in(a f), 
v(a) = 0, (3) 
p(B)v'(B) + bv(B) = 0, 


whose eigenvalues are given by the minimum principle 


f 
i (pb? + ad?) + 698) 


U = min =— y ——_—— 
f py? 
a 
the minimum being taken over all functions $ which satisfy the conditions 
(a) œ is continuous and piecewise smooth on [«, f], 


(b) (a) = p(B)$'(B) + boi) = 0, 


(c) ipii = [ ppn = = [oom = 0, where v,,...,U,-, are the eigen- 
functions corresponding to the first k — 1 eigenvalues. 

Any minimizing function of this class which is twice continuously differentiable on 

(æ, B) is an eigenfunction v, of the system (3). à 

The kth eigenfunction has exactly k — 1 zeros in (a. ff). 

Increasing p or q or b, or decreasing p. or taking a subinterval of [, f] increases the 

eigenvalues of (3). 


The eigenvalues of the system (3) separate those of the system (1) and vice versa. 


13.5 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 
The general solution of the equation is of the form 
u(x) = A cos 45x + B sin 23x, 
and the boundary conditions give 
A cos Ata + Bsin Ata = 0 (1) 
and 
A cos À* ff + B sin 4B = 0. (2) 


These equations have a nonzero solution provided that the determinant of the 
coefficients is zero, i.e. if 


cos Atay sin Jf — sin Ata cos AB = 0 
or 
sin AB — a) = 0. (3) 


Therefore the eigenvalues are given by 


à= hra = = n=1,2,.... (4) 

(f — a) 
The solution 4 = 0 of (3) is rejected since it gives A = 0 and therefore the trivial 
solution u = 0. The solutions corresponding ton = — 1, —2,... merely duplicate (4), 


and so are omitted. To find the corresponding eigenfunctions we solve (1) and (2) with 
4 = À, there will clearly be an infinite number of solutions, for which we have, in 
general, 

nme 
B-a 
from Equation (1), which is the same as — tan Af obtained from Equation (2). The 
eigenfunctions are therefore given by 


A 
g^ Ata = —tan 


w(x) = C sin ee cog EX +. cog TE sin I 
al? Er ET usa ean 
. nu(x — a) 
= Csin "X T 


where C is an arbitrary constant. The end point « has no special standing in (5); an 
alternative process would have given 


nnz(x — fj) 


u(x) = D sin ————., 


B-a 

which is the same as (5). 
Solution to SAQ 2 
The general solution is of the form 

A cos 43x + B sin Aix, 
and the boundary conditions give 

— AA! sin 2ta + Bj! cos Ata = 0 (1) 
and 


— Ast sin 148 + B4? cos 23 = 0. 
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If 4 and B are not both Zero, we must have 
A( — sin Ata cos 2+8 + sin 2f cos Ata) = 0, 
or 
Asin ANB — a) = 0. 
The eigenvalues are given by 


n?n? 


iie mer med s... 


We do not reject 4 = 0, since this does not yield the trivial solution in this case, In 
fact, corresponding to the eigenvalue 2) = 0, we have the solution 


Ug = A, (3) 
where A is any constant. 
For the other eigenfunctions we put 2 = 2, (n > 1) and solve (1) and (2), obtaining 
nn(x — a) 
p-a 
where C is arbitrary (or D cos [nz(x — B)/(B — a)] where D is arbitrary). 


To resolve the paradox of 4 — 0 being an eigenvalue, we follow through the argument 
of W: page 161 for our problem. Equation (36.4) is still true (2 being any eigenvalue 
and u a corresponding eigenfunction) with a suitable change in the limits, because 


2 
par | =0. 


x 


u,(x) = C cos 


The function q is nonnegative for our case—in fact q = 0, The difference lies in the 
fact that 


u(x)=0 a<x<f, 


is now a possibility; indeed u given by Equation (3) satisfies this condition. Therefore 
the numerator of Equation (36.4) can be zero and the eigenvalue 2 = 0 is correctly 
given. 


Solution to SAQ 3 


Suppose that there exist two linearly independent eigenfunctions corresponding to 
Ay, uj? and uf?. Let u = U be a particular solution of the equation 


(pu — qu + A,pu = 0, (1) 
for which U(x) = 1, but which is otherwise unrestricted. The existence of such a 
solution is guaranteed by the Existence Theorem (Unit M20/ 33, Existence and 
Uniqueness Theorem). Then, since ut!) and uf?! are by hypothesis linearly independent 


solutions of the second order differential equation, they form a basis for the solution 
space and there exist constants A and B such that 


U(x) = Aux) + Bul?(x) x € (a, ff]. 


But the right-hand side is zero when x = a, and the left-hand side has the value 1. 
Therefore we have a contradiction, and we conclude that any two eigenfunctions cor- 
responding to the same eigenvalue are linearly dependent. 
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Solution to SAQ 4 


Let 4, jt be two different eigenvalues, and u, v a pair of corresponding eigenfunctions. 
We carry out the procedure in W: pages 160 to 161; u and v satisfy 


(pu) — qu + Apu = 0, 
(pv) — qv + ppv = 0. 
We multiply the first equation by v and the second by u, subtract, and integrate over 


[x, f); after cancellation we have 


h p i 
Í [vo(pu'Y — u(pv')] + (A — n puv = 0. 


After integrating by parts (an alternative view to that of W: page 161, line 4) and 
rearranging, we have 


Li | 
A= n pub = LE - «)| . (1) 


We aim to show that the bracket is zero. Since A # 0 we can write 
1 
plou’ — ww) = gre + Bu) — u(Av' + Bv], 


which is zero at x = a, since u and v satisfy the boundary conditions. Similarly we 
may show that it is zero at x = fJ. Therefore (1) gives 


p 
(4— nf puv = 0, 


and since À # y, 
Li 
Í puv = 0. 
a 


Solution to SAQ 5 


Such periodic boundary conditions arise when, for example, 0 represents the angular 
coordinate in plane polar coordinates. When 0 = n we start “going round again”, and 
the boundary conditions state that the function and its first derived function (and, 
because of the differential equation, all its derived functions) are continuous across 
0 — x. The particular equation in our problem arises from separating Laplace's 
equation in polar coordinates (W: Section 24). 


The general solution is u(Ü) = A cos 230 + B sin 450, where A and Bare any constants. 
The boundary conditions give 


u(—n) — u(x) = 0 = —2B sin 4*2 (1) 
and ` 

u(—2) — u(n) = 0 = 242? sin Ata. (2) 
One solution is 23 = 0, which gives the eigenfunction 

u(0) = A, 
where A is any constant. 


If 2 # 0, then 4 = i? (n = 1, 2,3,...) satisfies (1) and (2). Both A and B may have any 
values, and the corresponding eigenfunctions are given by 


u(0) = A cos n0 + B sin nd. 


The simplest basis for the whole space of eigenfunctions corresponding to the eigen- 
value i? > 0 is 


{cos n0, sin n0}. 


Thus we see that the eigenspace has dimension 2. 
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Solution to SAQ 6 
The general solution is 
u(x) = A cos Atx + B sin Aix; 
the boundary conditions give 
u(0) 20-4 (1) 


and 
w(2n) — u(2n) = 0 = B(A* cos 232n — sin A32n). (2) 


paanan (2) gives 4t = tan(2n23), as required. There is an infinite number of positive 
igen alues, as well as 2 = 0, as can be seen by plotting a graph with axes 4, y, and 
sketching the line y = 4* and the curve y = tan(2n43). Do 


To find any negative eigenvalues we wri 
re write = =u (u > i = jut H 
for the roots becomes dts andà ius the equation 


j3 = tanh(2zj), 


since tan ic — i tanh c. The sketch Shows two values of y+ equal in magnitude and 
opposite in sign, giving a single negative root 2. 


A 


y 


Note that the existence of a negative eigenvalue does not contradict W: page 161, 
line — 11, since the boundary conditions in our problem are different and (36.4) does 


not hold. 


Solution to SAQ 7 


The first eigenvalue is 2, = n? and the corresponding eigenfunction n,(x) is A sin zx 
(A some constant). The Rayleigh quotient becomes 


1 1 
Í An? cos?ux dx — A?m? f 4(1 + cos2nx) dx 
0 Jo a 


1 ni 1 
f A? sin? nx dx af 4(1 — cos 27x) dx 
«0 o 


as required. 


Solution to SAQ 8 


Let i be an admissible function which minimizes the given quotient, and let yr be any 
admissible function. The minimum property implies that 


RIY] > ROW). 
where R is shorthand for the given Rayleigh-type quotient. We write the comparison 
function y in the form 

y - oc. E 
where $ is an arbitrary admissible function. and v is an arbitrary parameter. (Obviously. 
we do not in any way increase the degree of "arbitrariness" available by introducing v: 
we can in fact do without it, but the argument is simpler if we put it in. It is convenient 


to think of t being "small" in some sense.) Then. for any given $ the quotient R[y + th] 
has a minimum when t = 0, that is, 


d [ iad to? — (WU) + tht) 


dt 


li (w+ th) 


t=0 


f 2j'à' — Wal) t y? - war 2yo 
NEL st < D =0. (1) 


NF. Ln] 


Let the minimum value of the quotient be equal to jt i.e. 


1 
[ v? - war 
zo 


Ry) = ———7————- = as 


E 


0 
then from (1) we obtain 


[, we - nin = vae = 0 Q 


for every admissible $. We integrate (2) by parts, assuming that y is twice continuously 
differentiable: then 


0 


= f $^ + uj) + WDA) — WOPO — yo) 


1 
“0 


-f $^ + p) + WO) — WNO) — v'(0)9(). 


The terms at the end come to zero. since V'(1) — (1) 2 0 and $(0) = 0, by the 
boundary conditions. The rest of the argument is as in W: page 164. 


Solution to SAQ 9 


Suppose that 2,,..., 2,-, and t... t, have been found. Then to find A,, i, we 
proceed exactly as in W with p = 1, q = 0, p = 1; integrating over (— m, n) instead 
of (0. 1), we obtain, corresponding to Equation (36.9), 


[tw -wo =0, Aj 


where y is the minimizing function and $ an arbitrary admissible function subject to 
the boundary conditions of the present problem and to the k — 1 orthogonality 
conditions (if k = 1 there is no orthogonality condition to be satisfied) 


We integrate (1) by parts: for all admissible $ we have 


f [-v" — upe + [vol =0. 


=a 


13 SAQ9 10 
Now 


[vo] = WIR\b(n) — pinpin) = 0. 


because of the boundary conditions 


satisfied by @ a Bá Sim —: 
W: page 164. only the boundary Mae ie Oe a “= 


Conditions being different. 


We saw in SAQ 5 that the eigenvalues are 4 = (n= 0, 1,2. ) If n = 0 the ei 
space has dimension 1. and if n > 0 it has dimension 2, In the minimization piores, 
the successive minima delivered occur in equal pairs (after the first) The kai ei en- 
function delivered in any such pair is, of course, diflerent from the first. si (d 
required to be orthogonal to the first. oem? 


Solution to SAQ 10 


We proceed exactly as in W up to Equation (36.9). obtaining 


ii [owo + ghe — pid) = 0. 


Now we integrate by parts, so that 


1 
= I [VY — qi + np] + POYDA) — p()y'(0)9(0) = 0. u) 


and for every admissible $ we have $(0) = 0, so the last term vanishes. 


Now in order that (1) should vanish for every admissible @, it must vanish in particular 
for such $ as satisfy additionally the condition 


$0) = 0. 
Therefore 


$ 
f Bl(pw'y — qu + upy) =0 


for every ó for which $(0) = (1) = 0, and the argument in W: page 164, lines 5 to 10 
shows that y must satisfy 


QUY — qi + up — 0 in (0, 1) Q) 
with. of course. 
W(0) = 0. (3) 


Therefore, from Equation (1). we are left with 

PDPO) = 0 
for every admissible @. Since p(x) is positive for 0 < x < 1, this is possible only if. 
additionally. 

y'() = 0. d) 
Thus the minimizing function y is an eigenfunction of the system (2), (3). (4) with cor- 
responding eigenvalue jt. the minimum value of the quotient. 
For the second part, we consider the system 

V" + 2p =0 in (0, 1), 

w(0) = V1) = 0. 
For this system, the result just proved shows that for any smooth function $ satisfying 
Q0) = 0. 


where Z, is the smallest eigenvalue. It is easy to confirm that the smallest eigenvalue 
hy E g i : 
is 1x? (see. W: page 68). and the result follows. 
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Solution to SAQ 11 


We can rewrite the differential equation in two ways, 


2 


ea) + dru =0 (d) 
dr? 

or 
ife z) + Aru =0.- o (2) 


(a) The first way enables us to get all the. solutions: these are given by 


cos.A*r sin 2ir 
u(r) =A + . 


r r 


The condition "u(0) is finite" gives A = 0. The condition "u(a) = 0" gives the 
eigenvalues 
24,2 
AE n=1,2,.... 
di 


The solutions (eigenfunctions) are therefore given by 


ur) = uui n=1,2,.... 
r 


(b) The second form, Equation (2), is the self-adjoint form of the given equation, 
and leads to the minimization problem: 


minimize 


f r (ert)? dr 


A ee 
f Pte o 
0 


over all continuous and piecewise continuously differentiable functions $ on 
(0, a) (in particular, finite at r = 0), satisfying (a) = 0. 


Solution to SAQ 12 


Elh: d) 
a) R = . 
(a) R[9] Hib.) 
(b) We confirm the inner product axioms one by one, 
H is 
symmetric: Clearly H (v, w) = H(w, v). 
1 
linear: H(u, ab + Bw) = Í pu(a» + Bw) 
o 


= aH (u, v) + BH (u, w) 
for all a, BER. 
positive-definite: H(u,u) > 0 unless u = 0 (since p > 0), 
and 
H(u,u) = 0 when u = 0. 


(c) Eis 
symmetric: Clearly E(v, w) = E(w, v). 
1 
linear: E(u av + Bw) = f [pu'(av' + fw) + qu(av + fhw)] 
o 


= «E(u, v) + BE(u, w). 
1 
positive-definite: E(u, u) = Í [pu + qu^), which is greater than zero for every 
o 
nonzero u provided that q > 0. Also when u = 0, E(u, u) = 0. 
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SAQ 12-13 


(d) Fora minimum, there exists an admissible function w such that 


EW + zx + vy). EQ.) 
HUE und) Hypa e SY É 


for all real rand for all X Tor which (0) 
the expression on the left ism 
term (using the symmetr 


whi (0) = x(1) = 0. Given any such function pa 
attains, its minimum when t = 0. The expansion of this 
y and bilinearity of E and H)is 


El V) + EN. 7) +t Ely, y) 
Hb, V) + 2xHJ: 7) C EQ, yy 
and the derivative at. = 0 is‘found to be 
HW, W2EW. y) — EW. W2H WW, 2) 
UP 077 


This is zero when 

E(W. x) — uH (p, xy = 0, 
using (1), and this must be true for every x. Substitution of the integral repre- 
sentations of E and H gives Equation (36.9) in W again, and hence Equation (36.10). 


(e) Let dp bea solution to Problem B. We wish to prove that it must also be a solution 


of Problem A. We make the hypothesis that dy is not a solution of Problem A, 
and show that this leads to a contradiction. 


Let $, be any solution of Problem A. Then it is scen easily that any multiple of 
$4 is also a solution of Problem A. We now define 


Da = aha, 


where a = [H (ha, 64))7 4; ®, is therefore also a solution of Problem A, with the 
property 


H(9,.0,) = PH lha, 0) = 1; (2) 
€, is therefore an admissible function for Problem B. Now, 


Elpo. dp) SN E(®,, $4) 
H (hg, o) | H(b,. D4) 


since $, is, by hypothesis, not a solution of Problem A. But 
H(p. dg) = 1 

by definition, and in view of Equation (2), (3) becomes 
Elp. Pp) > Ela. Da). 


This contradicts the fact that y is a solution to Problem B, and therefore our p 
must be a solution to Problem A. 


Solution to SAQ 13 
We require the eigenvalues and eigenfunctions of the system 
u” -iu-0 in (0, 1), 
u(0) = u(1) = 0. 
2... ) and u,(x) = sin nax. Then, by (36.14). 


2. sin mx sin nny 
G(x.) = — — 


1 
! we? [ sin? mé dé 
0 


2 £ sin nax sin nny 


x 
EET E 
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Solution to SAQ 14 

The functions given are the eigenfunctions of the system 
u.-ii-0 in (0. 1). 
u(Q) = u(1) = 0. 


which is of the type considered in this section. 


Solution to SAQ 15 o 
We have 
Ak € Aye 


and uj, ,. u, have respectively exactly k and k — 1 zeros in the open interval (0. 1). 
and are also zero at the end points. By the Separation Theorem (with 7 = Z,. 
A= AGaGaUC = Up. U m p4 1). there lies at least one zero of t, , on the open interval 
between each two successive zeros of up. Since there are k such intervals, the k zeros 
of u+, on (0, 1) must be those on the open intervals, exactly one on each interval. 


Solution to SAQ 16 


Apply the Separation Theorem with 4 = 4 and (a) u = u* v = u** (b) u — à 
The result follows since u* and u** are linearly independent. 


Solution to SAQ 17 


The eigenfunctions are displayed on W: page 170, line 18. The zeros x,,, say, of X, 
(n = 1,2,...) occur when 


log(1 + Xn) 
log 2 


g2 


= Mt 


m being any integer. and 0 < x,, < l: ie. when 


for some integer m. 
Now. since n > 0, 

m 

27 -1>0<em>0 
and 

m 

2n- |<lem<n. 


Therefore 0 < x, < 1&0 < m « n. so there are exactly n — 1 zeros, given by 
m= 1.2.0... n=l. 


Hi. 
For the second part the zeros of X,,,., are given by 2^*! — I {p = 1.2,..., 5). There- 
fore we need to show that 
al i 2 2 " 
O0c«2*1-1«2»—1«2*i—-|1zc«2:—1«..«25i—1zl, 


which is obvious. 
Solution to SAQ 18 
Consider the boundary value problems 
(pu) — qu + ¿pu — 0 in (0, 1). 
u(0) = u(1) = 0. (1) 


p(x) = 1 + x7. g(x) = x, p(x) = 1 + x2, 
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PDE 13 SAQ I8 19 
and 
(quyY-dqi-2zpi-0  in(0. 1), 
(0) = #1) = 0. 
Firstly, let 


(2) 


P=1.9=0 ps2: 


ji >A ü 5 
E 2 p. q > ij and p < p, and the two problems are not identical. By the 
onotonicity Theorem we find that the eigenvalues of the two problems are related by 


AX A. (3) 
Problem (2) becomes 
i” + 277 =0 in (0. 1), 
(0) = i(1) = 0, 
with eigenvalues 
Ay = dee. 
Therefore, by (3), 
ij eA. (4) 
Similarly, let us compare with (1) the problem (2) with pe 2, ĝ = 1, p = l; then 
20" — ū +77 =0 in (0, 1), 
or 
i’ +3- Ni =0  in(01) 
and 
ü(0) = i(1) = 0. 
The eigenvalues are given by 
i8, — 1) = Ez 
or 
A, = 2k? n? + 1. 
By the Monotonicity Theorem. 2k?z? + 1 > A,. Therefore, from (4). 


then? eA < 2k?n? + 1. 


Solution to SAQ 19 
(a) We proceed as in W: page 172. Choose the trial function 


0 aix<g, 
P(x) = 4B x) ZS NSB. 
up B<x<Bp, 


where t, is the first eigenfunction of 

(pi) — qt + up = 0 in (3. ff), 

EE) = PEB) + bit) = 0. 
Clearly. @ is continuous and piecewise smooth on (x. fj]. We find that provided 
q = 0, 


RE 


D - 
[ lpo? + q4?] + bio? f [PE + qe] + ble AYP 


—— —— = fi,, 


since olf) = v). p > p. d 2 q. b = b and p < p. The proof then continues 


as in the text. 
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PDE 13 SAQ 19-20 


(b) Compare the given system with the system 
v" + 42v=0 in (0, 1), 
v(0) = v'(1) = 0, 
where we have decreased p and increased p. The lowest eigenvalue is given by 


24, = (n 


or " 
My = ar. 

Therefore 
A, > dn’, 


by part (a) of the question. 


Solution to SAQ 20 
We proceed as in W: pages 173 and 174. Suppose that u, has / zeros, at « = Xo, 
Xp Xi, f. For m 2 1,2,...,1H let 
[63 Xar RRS X 
ye) = { K(x) 1 
0 elsewhere, 


where x, = f, and let 

$(x) = p Cy uf), (1) 
where €... C are any constants. Then, as in note (ii) of Section 13.3.2, 

[in^ + 09% + ito 


Xm L 
&[^ tm? + aud) e Y, atop 


Xm m=1 


i 
X 
mzl 
I Xm Nm 
PE f su] 7 [ [pii = oul} + beth)? 


Xm 


i 
PEUL) — Y. cn gl (puny! — qu] + bet pu)? 


m=1 M 


(noting that, this time, u,(fj) # 0) 


Xm 


l 
= eu) pui) + bu) — Y. cs ullpu) — quj). 


m=1 Xm-t 


The first term is zero, by the boundary condition at fj; therefore, using the differential 
equation, this expression becomes 


l 
ag 2. Cn 


m=1 Xmas 


Xm 


. rj 
pui = af po. 
Hence, 


B 
[ in? + 0071 + broo 


r pp? 


=A, (2) 
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20-21 


where ¢ is defined as in (1), and we have not so fa 
any way. The | constants {c,,} m 
equations 


r restricted the coefficients {¢,,} in 
ay clearly be chosen to satisfy the | — 1 linear 


n 
Í pou, = 0 r=1,2,...,1-1; (3) 


that is, we can arrange for $ to be orthogonal to the first | — 1 eigenfunctions of the 
problem, without disturbing (2). 


Now we know that if we minimize the quotient in (2) over all admissible functions $ 


which also satisfy Equations (3). we obtain a minimum equal to the /th eigenvalue 2,. 
Therefore 


2 S ds 


ISk, 
w r the number of zeros of u, on (a, B), which is | — 1, must be not greater than 
Solution to SAQ 21 
Let v be any solution of the equation 
(pv') — qv + ppv = 0, 
with 
vla) = 0 
and let u satisfy 
(pu — qu + Apu = 0, 
with 
u(a) = u(B) = 0. 
Let | be an integer >1, and suppose that 
Any SH <A. (1) 
We shall show that v has exactly | — 1 zeros on (a, fj). 


By the Separation Theorem, since 4,_, < u, there is at least one zero of v between 
every consecutive pair of zeros of u,_,. But by the Oscillation Theorem, there are 
exactly | zeros of uj. , on a < x < fj, the end points being zeros. Therefore, if m is the 
number of zeros of v on (z, fi), 


mzalc-1. (2) 


Again, by the Separation Theorem, since jt < ,. there is at least one zero of i, between 
each consecutive pair of zeros of v: also, there are exactly ! — 1 zeros of u on (a. fj). 
Therefore, since v(x) = 0. 


mt&l- Il. (3) 
(2) and (3) together give 

i-1&mst&l-1, 
that is to say. 

mz-l-— I. (4) 


Now let p, be an eigenvalue of the system (38.6). It cannot be equal to A for any L 
since in that case n, would be an eigenfunction of both problems and satisfy both 


u,(p) = 0 
and 

pesto) + be.) = 0. 
That is to say. 


eff) = nf = 0. 
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But by the Existence and Uniqueness Theorem this implies that v, = 0. Therefore ju 
must satisfy 


Apa € d € A 


for some l, It follows from (1) and (4). and the fact that v, has exactly k — 1 zeros on 
(x. ff). that | — k. 


Solution to SAQ 22 


We compare the systems 


D" — quU guo = 0, r(0) = e(1) 2 0: (1) 
u" — qu + Au = Q. u(0) = u(1) = 0: 
= Goin + vw m0, — w(0) = wl) = 0: (2) 


where gmax: Gmin are the maximum and minimum values of q on (0, 1]. 


The Monotonicity Theorem shows that, for each k, 


ve S A, S pt. (3) 


with equality if and only if q is a constant function. But the eigenvalues of (1) are 
given by 


12-2 
Hy = Qmax = CR. 


and those of (2) by 


j pinza 
Wk 7 min = KET? 


Therefore, by (3), 


or 


Fn + din S Ay SAP quas 


n2 nin ^k ^k < <2 + max, 


ur ca) k? 


Now take the limit as k ~ %. 


Solution to SAQ 23 


(a) The eigenvalues of the systems with Pmax or Pmin in place of p are given by 


(b) 


40 


k?n? k?n? 


B= apu. (B 0)? pus 


Since 


Pmin S P(X) S Pmax 
the result follows from the Monotonicity Theorem. 
The kth eigenfunction u, has exactly k — 1 zeros besides the two end points. On 
each interval (x;, x; 1). uy must be the first eigenfunction for this interval, with 


the corresponding eigenvalue equal to 4,. The result (a) applied ro this interval 
therefore gives 


2 2 
Ds " n? 

E H S ALS — Ly . 

(Xizi — Xi) Pimax (Xiri — Xi) Puis 


It follows (after taking the square root) that 


x 
p S pia Xis1 — X). 
^k 


Phmin(Xit1 — x) € 


fori 20,1,...,k— 1. 


PDE 13 SAQ 23 
(c) Sum the inequalities obtained in (b), and put 
Nit — Xp = 6;, 
so that 


k-1 kn ket 
> Phminds S -F S > [NC 
i-0 k 


Now assume (we have not proved this) that 


lim My, = 0. 


o that all the intervals (Ni, Xj, 1) tend to zero in length. Then we obtain in the 
imit 


4 
o kn 
Í p(x) dx = jim n 


by the definition of the integral. This may be rewritten, after squaring, as 


_ A n E 
lim @ = EM qu (1) 
fi P(x) «| 


which gives an estimate for the eigenvalues for large k. 


The final result indicates that, for large k, the eigenvalues tend to those of a certain 
simple harmonic problem 


u” + Apu = 0, ulo) = u(B) = 0, 


n 2 
p= jii ZI i 


this is a kind of average value of p on the interval, though not, perhaps, the expected 
one. 


where 


As an example, consider the variable string problem of W: Section 37, with 


k?n? 1 
p(x) = T+? and A, = (log 27 + zi 
Equation (1) predicts 
Ax n? n? 


eC d f = (log 2)?” 
is (1 + x) | 


which is clearly correct. k does not have to be very large for this estimate of the eigen- 
values to be quite good: for example, the error for k = 2 is 0.3*;. i 
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PDE 14.0 
14.0 INTRODUCTION 


This unit is concerned with Bessel functions and some typical boundary value prob- 


lems in which they arise. Besse] functions arise as solutions of a certain second-order 
ordinary differential equation, known as Bessel's equation. The interest in such func- 
tions in a course on partial differential equations is due to the fact that Bessel's 
equation appears, for example, when we apply the separation of variables technique 


to time-dependent problems such as the two-dimensional wave equation and heat 
conduction equation in polar coordinates. 


Bessel functions belong to a group of functions which go under the general name of 


special functions, If you want to get some idea of the extent of the subject of special 


functions, borrow a copy of Handbook of Mathematical Functions edited by M. Abra- 
mowitz and Irene A. Stegun. 


Following the treatment in W, we look firstly at some properties of Bessel functions 
and then consider, by the use of illustrations, just how they arise in physical problems, 


PDE 14.1.1 


14.1 BESSEL FUNCTIONS OF THE FIRST KIND 


14.1.1 Solution in Series 


Notes 


(i) W: page 179, lines 2 to 7 
W: Section 39 is not included in the course and consequently this reading passage 
is your first contact with this particular boundary value problem. Completeness 
of the eigenfunctions will be assumed. 


(ii) W: page 179, Equation (40.3) 
Bessel's soris (40.3) is often written as 


du 
EL es "hu d — mju - 0. 


We assume that m is nonnegative. 


(iii) W: page 179, line 13 
The problem is to find constants o, dy, a}, à5.... Which ensure that the series 
satisfies the differential equation. The condition ay 0 is not a restriction but 
merely a way of stating that « is the first index in the series. Since the equation 
is of the second order we expect two such series solutions. A series of the form 


æ% 
Y by" 
k=0 

is called a power series 


(iv) W: page 179, lines 14 to 16 
This note fills in some gaps in the manipulation. Term-by-term differentiation 
of the series gives 
d © 
PIU = P (a + n)a,r**"7! 
so that 


os 


= 2 a + n)a,t**". 


He 
Differentiating both sides again with respect to t, we find that 

d du o 2 

iito m (a + natt, 


We now substitute the series for u and d/dt(t du/dt) in Equation (40.3) and obtain 


E E E 
e^ Y @+ nat" — mY a+ E sn =0. 
n=0 n=0 n=0 


The last series on the left-hand side may be rewritten as 


and line 16 follows. 


(v) W: page 179, lines- 11 to —7 
The choice 2^ "/m! for ay is conventional and the notation J,(r) is introduced 
with this value of ay. Note that all the remaining coefficients are completely 
determined. Since a, = 0 all coefficients with odd suffixes must vanish by the 
recursion formula 


m" dn-2 Loo 45-2 
(m + n — m? (2m + n)n 


PDE 14.1.1 
For the other coefficients we have successively 


ao 27m-2 


2(m+1).2 “ime hr 
a gms 
a4 PN 
2(m + 2).4 ^ 2(n 4 2)! 
and, by induction, 
» 4-2 = (= 127-2 
2(m + k).2k — k(m- E 


a, = 


Equation (40.4) follows by substituting these coefficients into the series. 


Jmis bounded at t =0 because m is assumed to be nonnegative. It is of course the 
boundedness condition which forbids us the choice « = —m. 


(vi) W: page 179, footnote " 
If m is not an integer, the convention for ag becomes 
27^ 
a = ——— ——. 
9 "Tim + 1) 
where the Gamma function F is defined, for m > 0, by 
EI 
T(m) = f e7*x"7! dx, 
0 


When m is a positive integer, T(m + 1) = m! (See SAQ 3.) 


(vii) W: page 179, line —5 

The formal construction of a series solution of a differential equation does not 
ensure that the series converges for all (or any) values of t. Before we can really 
define J,,(t) by Equation (40.4) we must examine for what values of t the series 


converges, The ratio test is a test for convergence. It states that if, in the series 
x 


Y Cys Ck x 0 and 
=0 


k 
lim |-.| <1 
kee] Cy 
E 
then the series $, |c,| converges. Now a series J` cp for which 5 |c,| converges 
k=0 


is said to be absolutely convergent, and there is a result which states than an 
absolutely convergent series is convergent.* 


In the series (40.4) put 
(- 1 (1o *2* 
Lm kl(m + K)!, ` 
whence, for t # 0, 


; 
- lim ai 


Ciel - 
ko (k + D)(m +k + 1) 


C 


lim 


kx 


for all t + 0. Thus, using the ratio test. and noting that the series obviously con- 
verges for t = 0, we see that the series converges for all t. For nonnegative values 
of m the Bessel function J,, is defined by the power series (40.4). 


General Comment 


You may wonder what happens if we consider the unbounded solution of Bessel’s 
equation. Put x = —m and carry the steps through formally. Then apparently 


< (gt 
egit à KIT(=m +k + 1l) 


*These results may be found in M. Spivak, Calculus (Benjamin. 1973). pp. 393 and 397. 


for any nonnegative m. But the terms of the power series are not all defined; for 


example, for k = 0 and m = 3, 


T(-m +k + 1) =F (-3) 


and the Gamma function cannot be defined by the integral formula in note (vi) since 
the integral diverges. However, provided 0 « m < 1 we get a satisfactory definition 
for J_,,. (In fact there is a different definition of the Gamma function which applies 
in the extended domain containing negative arguments, although we do not intend to 
go into this here.) What we do is look into the case « = —m by using the recursion 
formula 


(n — 2m)na, = —a,-5. 


Since a, = 0 again all the odd coefficients vanish. If m is not an integer then we can 
make some choice of ay and obtain another solution, independent of J,,, which will 
be unbounded as t ~ 0. If m is an integer, it follows from the recursion formula that 
055,-2 = 0, and by successive application of the formula we conclude that 


a3, =0 r=0,1,...,.m— 1. 


Thus, if m is an integer, we do not obtain a further independent solution. (Although 
ao = 0 is incompatible with the defining condition ay # 0 we could investigate this 
case further. The first nonzero coefficient would then be a;,. Assign the value 
(—1)"27"/m! to azm- It then follows that we could define 


co piggt 


Jon) = 3, ie n 

A change of the dummy variable from k to j by k = j + m now gives 
eo (—1y* "(Ay *2j 
-w(t) = ———————— = (-1)" 

Foul) = Y MW = O D, 
confirming that no new solution has been obtained.) 
SAQ 1 
Show that 

dJo(t) A 


d 7-0. 


(Solution on p. 20.) 


SAQ 2 
W: page 181, Exercise 2 
(Solution on p. 20.) 


SAQ 3 
(a) Using the definition of the Gamma function show that 
T()21, TQ)=1, r(3)=2. 
(b) Show that 
T(m + 1) = m (m) 
for all m > 0. Deduce that if m is a positive integer T(m + ])2m!. 


(Solution on p. 20.) 


PDE 14.1.1/14.1.2 


SAQ 4 


Given the result 


co 

- 14 
Í e dt = ixi, 
0 


deduce that I (4) = zi. 


(Solution on p. 20.) 


SAQ 5 


Supply the details to establish Equations (40.5) and (40.6) in W: pages 179 and 180. 


(Solution on p. 21.) 


SAQ 6 
Verify that 


ten 
Jo(t) = Hi cos(t sin 0) d0. 


(Solution on p. 21.) 


14.1.2 Zeros of Bessel Functions 


Notes 


(i) 


(ii) 


(iii) 


(iv) 


W: page 180, line 11 
That 4," increases with m follows from the Monotonicity Theorem on- W: 
page 174. 


W: page 180, lines 13 and 14 

Since J,, satisfies Bessel's equation it is certainly continuous for t > 0. Therefore, 
between any two successive zeros the function t^"J,, must have at least one 
stationary value. At this point d[t^"J,,]/dt must vanish; thus by (40.5), Jm+1 
must have at least one zero between successive positive zeros of J,,. We prove, 
by contradiction, that J,,,, has just one zero in this interval. Suppose J,,, , had 
more than one zero in (j,", j",). Then 


d 
diu aa] = 0 


between two such zeros, and from (40.6) with m replaced by m + 1 it would follow 
that J,(t) = 0 for some t in (j,?, j£). This would contradict the assumption that 
the original zeros of J,, which we chose were successive. Thus the zeros of In and 
Jm+1 are interlaced. 
W: page 180, Equation (40.9) 
This problem is not equivalent to the original eigenvalue problem (40.1) and (40.2). 
Solutions for which w(0) = 0 do not necessarily satisfy the conditions 

u bounded 

lim xu’ = 0 as x ~ 0. 
However, we now know the properties of Jat AX) as x ~Q. 


W: page 181, lines 4 to 6 
The notation 


is often used to mean that 


In this case we often say that g(x) is an asymptotic form for f(x). 


(v) W: page 181. lines 8 and 9 
In the terminology of IW: page /62. line —/ the eigenfunctions are 


(x) = Ju Lum) 


and the weight function is p(x) = x. The integral 


1 
f XI phx) dx 
Jo 


appears in the denominator of the Fourier coefficient c,. 


(vi) W: page 181, line 11 
In the right-hand side of the equation J2(t) means [J,()]*. . 


The following example shows you how to construct the Fourier-Bessel series (as it is 
known) for a given function. 
Example 
Find the expansion of | in the interval (0, 1) in terms of the eigenfunctions Jo / A, 9 x). 
We require the coefficients A, in the series 

eo ——— 

1-2 Y Apdo fA) — xe(0 1). 

ket 
We multiply both sides by xJo(,/2,9)x) and integrate from 0 to 1. By the ortho- 
gonality property of the eigenfunctions we are left with one term on the right-hand 
side. Thus 

i = 1 = 
f xA Lx) dx = A, f xU LL Mx}? dx 
Jo 
= FAL, 


by W: page 181, line 17. The integral on the left becomes 


1 -- po pva 
i XJ, x) dx = = f tJolt) dt 
0 4 o 


1 ain 
- gs [n] by Equation (40.6) 


Hence 


M 


qim oye 


and 


SAQ 7 

W: page 181, Exercise 3. 
(Solution on p. 22.) 
SAQ 8 

W: page 181, Exercise 4. 
(Solution on p. 22.) 
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Extensive tables of Bessel functions for diffe 
following sketch shows graphs of Jolt), F(t) 
their oscillatory char: 


rent orders exist (see Bibliography). The 
1 and J,(t) in the interval 0 < t < 20, and 
acter (cf. the Oscillation Theorem) is clearly visible. 


—0.5F 


The table below gives the first few zeros of Jo(t), J,(t) and J,(1). Here j,"" represents 
the kth positive zero of J,(t). 


k Ji i? Je 

l 2.405 3.832 5.136 
2 5.520 7.016 8.417 
3 8.654 10.173 11.620 
4 11.792 13.324 14.796 
5 14.931 16.471 17.960 


We shall see in Section 14.2 that the Bessel functions of integral order arise as solutions 
of an equation which is derived by modelling oscillatory systems in two dimensions. 
This is the physical background to the oscillatory nature of the J,,,. You might wonder 
whether these can be related in a simple way to the more elementary trigonometric 
functions which have the same oscillatory property. In general this is not the case: 
there is however one special set of Bessel functions for which there is a simple relation- 
ship with the trigonometric functions. This is the set of Bessel functions of half-integral 
order, i.e. 


Just m=0,1,2,.... 


(The functions 


4 
2M 
i Ja s a) 


j| 
are known as spherical Bessel functions; they arise in the solution to the wave equation 
in spherical polar coordinates. We shall not however give a derivation here.) 
We now obtain the relationships by deriving firstly the one for J, . and then using W: 


page 179. Equation (40.5) to determine the others. J, is the solution of 


u 
'y-—-tu-0 
(ru) at u 


which is bounded 4t 1 = 0. 
Make the substitution 
w(t) = uni: 
then, as in Equation (40.9), w satisfies 
w'+w=0, 
w(0) = 0. 


This has the independent solutions sin t and cos t and, since w(0) = 0, we must discard 
cost and get w(t) = A sint as the solution. Thus we can argue that 
Asint 

HO 


JÌ = wo = 


We can settle the value of A by taking the limit as t ~ 0: 


lim ap = xm from Equation (40.4), suitably adjusted 
whilst 
int 
lim 225, 
090 t 


Matching these leads to 


1 
21T (3) 
since T(3) = 3T (3) and T (4) = J/z (see SAQ 3 and SAQ 4). Thus, finally, 


à 
JA) = B sin t. 


24 


A 
x 


[You may also derive this result by writing out and reorganizing the-series expansion 
W: Equation (40.4), and recognizing the Taylor expansion of the sine function.] 


To determine the formula for the other Bessel functions of half-integral order we 
write W: Equation (40.5) as 


1 ld[i 
pe Init) = E Jino] 


We may use this formula recursively to obtain, by induction, 


1 1dVTl1 
— 4, 2[--—]|— $ 
ghi) = [75 s) [m 20] 
In particular, for m = 4 we have 
1 ldYV[1 
phat) = -74] an0] 


or 


20 "Isi 
Fray) = "E [- ii (= t 
E un t 


We have indicated in Section 14.1.1 that for m » 0 we could obtain another (inde- 
pendent) solution of Bessel's equation by the power series method, choosing a= —m. 


Provided m is not an integer we obtain a new solution, independent of J,,, which is 
unbounded at the origin. For m = 4 the series 

S C DG 

Jat) = Y Trea 
ceo KIP(k + 3) 


is obtained. For m — n + $(n = 1,2,...) the expression l'( m + 1) is not defined. 
We can however define J_,, recursively using Equation (40.6). It is possible to show 
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that the function J_,, so constructed isa solution of Bessel's equation of order m and 
is not linearly dependent on J, 


m: In the next SAQ we ask you to derive an expression 
for J „-4(t) when n = 0, 1.2, 


SAQ 9 
(a) Show that 


+ 
J-i) = B cost. 


HINT: J., and J, are related by W: page 180. Equation (40.6). 
(b) Deduce that 
213 [1 d [cost 
J-,-40) = (=) [= —| |——— 
a e| | a | t 


T 


for n =0,1,2,.... 
(Solution on p. 23.) 
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14. TIME-DEPENDENT PROBLEMS IN 
TWO DIMENSIONS 


14.2.1 Vibration of a Circular Membrane 


Note 


(i) W: page 182 
This is the membrane problem which first appeared in W: page 48. The mem- 
brane, which is stretched across a fixed plane circular wire, is given an initial 
transverse displacement u(r, 0,0) = f(r, 0) and released from rest. Subsequently, 
the membrane vibrates freely. 


Q 


Show that the solution to the vibrating circular membrane problem (W: page 182) 
has the form (41.3). 


(Solution on p. 24.) 


SAQ 11 
W: page 185, Exercise 1(b) 


(Solution on p. 24.) 


SAQ 12 


A membrane is stretched across a circular frame of unit radius. The transverse dis- 
placement u(r, t) satisfies. 


Ou 

2x72. 
ay = VU. 
or 


The membrane is released from rest in the position u(r, 0) = «(1 — r°) where s is 
small and r is the distance from the centre. Find the subsequent displacement of the 
membrane. 


(You will need the result of SAQ 7.) 


(Solution on p. 25.) 


14.2.2 Sound Waves in a Tube 


We begin with a generalization of the discussion in Unit 1, The Wave Equation in 
which the one-dimensional wave equation was derived for the propagation of sound 
in a gas. We show how the wave equation in two or three dimensions follows from the 
equations of motion, which are: 

Equation of Continuity 


op i _ 
E + div(pv) = 0, (1) 


which expresses conservation of mass; (See Appendix.) 


Momentum Equation 


= —grad p, (2) 


Ov s 
ae + (v-grad)y 
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assuming the gas is Subject to no external forces [this equation may be interpreted 
componentwise. i.e., for any fixed unit vector e 


ev, 


PW [a erani] = —e-gradp* 


where v, = e+ v is the component of v in the direction of e]; (See Appendix.) 


Equation of State 
=p", 
p p (3) 


where a and y are constants depending on the particular gas. (See Unit 1, Section 1.1.3.) 


These equations provide relations between the 


ressure p(r, t), the densi " 
the gas velocity vr, t). P PA UR density AC RUM 


Suppose now that a sound wave passes through a gas at rest with initial pressure 


P(r, O) = py and density p(r,0) = Po. Equation (3) must be satisfied in the initial 
state so that 


Pat pP = poll + s), where we assume as part of the approximation that s is small. Thus, 
rom (3) 


p ; 
pal + s)” = poll + ys), 


p= 

p 
where we have retained just the first two terms in the binomial expansion of the right- 
hand side. The basis of the approximation is that s and v and their derivatives are 


small so that second degree terms containing them may be neglected. Thus Equa- 
tions (1) and (2) become 


ô i 
x t div v — 0, (4) 
ðv 
Pog = —grad p ~ — poy grad s: (5) 
Assuming that p (and hence s) is twice continuously differentiable we have 
ov as 
9v _ _ a's 4 
div ET ad from (4) 
= — P diy grads. from (5) 
Po 
We set c? = ypo/po and obtain 
1 @s 
*s =z =0 
Vis ~ agp 7% 


i.e. the density satisfies the two- or three-dimensional wave equation. 


Since the gas starts from rest it may be shown that, as a result of Equation (5), there isa 
scalar field ® such that 

v = —grad®. 
Note that ® is undetermined to the extent of an arbitrary function of time. We call «b 
the velocity potential. 
Rearranging (5) we find that 

graa [s = z2) =0, (6) 


Po et 


whence 


An arbitrary function of time could appear from (6) but we absorb that into &O/t. 
Finally substituting the expressions for s and v into Equation (4) we get the wave 
equation 

e 


L = div grad b = V^. 
rai" 


The number c is a constant which measures the rate of propagation of disturbances 
and for this reason is known as the speed of sound in the gas. 


Suppose air is enclosed by a long hollow tube of circular cross-section with radius a. 
The tube as a whole vibrates rigidly perpendicular to its axis with velocity ugcos at, 
that is, it oscillates harmonically. Let us find out how the air inside the tube vibrates. 


Up COS at 


Fora long tube it is arguable that, away from the ends of the tube, the velocity potential 
is independent of distance along the tube and varies only with time t and position 
relative to the axis of the tube and the direction of vibration. Polar coordinates are 
suitable for this problem because the boundary is a circle in the cross-section. Thus 
«p satisfies 

Op 10> 100 1 


= t- tha (Kr K<a0<0<2n 
e ro ra co? i ! 


Plho=0 = loca. 


and $ is bounded as r ~ 0. How does the air immediately adjacent to the wall behave? 
We take no account of viscosity; so we suppose that the air may slip freely parallel 
to the wall of the tube, but must have the same velocity normal to the wall as the wall 
itself. Let 0 = 0 be the direction of vibration. The normal component of the wall's 
velocity at points with angular coordinate 0 is ugcos ot cos 0, whilst the normal 
velocity of the air is — 60/dr. Thus 


[i 
-—— = ugCOS at cos 0, 


The boundary conditions suggest we try a solution of the form 
D = f'(r)cos 0 cos at 


where, after substitution into the wave equation, f satisfies 


1 
J" + Fu -3/ü-EHft)20 O<r<a, 


S'a) 2 — ug. 


and k = a/c. This is Bessel's equation of order 1. The solution which is bounded at 
r = Ois 


fir) = AJ, (kr). 
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The constant 4 is determined by the boundary condition 


kAJ\(ka) =-ug > 45 —-0_. 
! * KJ' (er) 
The required solution is therefore 
tot J (Kr) 
@= — 2 ENM 
v Jika) cos ( cos ct. 


Since J, is an oscillating function it follows that Jı also oscillates about zero; in other 
words J, has an unbounded number of zeros. If ka coincides with any one of these 
zeros then the solution breaks down. In the neighbourhood of these zeros the ampli- 
tude of the wave is very large: the linear theory of sound waves is then no longer 


appropriate and the assumptions made in the derivation of the wave equation must 
be re-examined. 


14.2.3 Heat Conduction in a Bar 


We complete this section by looking at an example in heat conduction. 


Consider a long bar of circular cross-section and unit radius which has an initial 
temperature distribution which is radially symmetric and does not vary along the bar. 
Suppose the boundary temperature is maintained at a constant value, say zero. The 
bar is supposed to be long so that we can ignore end effects. What is the temperature 
distribution at time t over the cross-section of the bar? Since the temperature depends 
only on the radius r and time t, the heat conduction equation becomes 


ĉu ĉ 
e , a i ] 
where u(r, t) represents the temperature. The boundary and initial conditions are 


respectively 
u(1, 1) 2 0 t>0 


O<r<lsr>d, 


and 
uO =f) O<r<] 


where f is a prescribed function. We also require that lim u(r, t) = 0 and that u(r, r) 
tn 
be bounded as r ~ 0. 


To apply the separation of variables method we look for solutions of the form 
R(r)T(r). Thus 


] df dR _ lt dT_ 
rR de dr kT dt o 


say. The functions R and T satisfy the ordinary differential equations 


ORY + ArR =0 Q«r«l. R(0) bounded. R(l)- 0, 


and 
T' +A4kT=0 t>0, lim T(t) = 0. 


(mn 
The second problem has the solution 
Tit) = e7. 


We require the solution to decrease with time and so we impose the ont 40. 
The first equation is Bessel's equation of order zero with a solution Jo / 2r). The 
other independent solution of Bessel's equation has not been investigated when the 
order is an integer but it has a singularity at r = 0, whereas Jy is bounded at r = 0. 
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The boundary condition R(1) = 0 is satisfied if we choose A so that 
Jo /4) = 0. 


The solutions of this equation are written conventionally as 4, (n = 1,2,...). Thus 
the full solution is 


© 
u(r.t) = Y, Ago / Annem an, 
n= 


The constants A, are determined by the initial condition u(r, 0) = f(r), i.e. 


E 
fe) = Y AJ G/ Ai). 
n=l 
Hence A, are the coefficients in the Fourier-Bessel series for f(r). 


SAQ 13 


Solve the heat conduction problem in a cylindrical bar of unit radius with the initial 
temperature distribution given by 


Sf =1— 9. 
(You will need the result of SAQ 7.) 


(Solution on p. 25.) 
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14.53 SUMMARY 


This unit introduces the Bessel function of the first kind Jm and some of its properties. 
The Bessel function is defined as (a power series which is) the bounded solution of 
Bessel's equation of order m 


2 
aUu du 


xL e opu TO EE LT Ce 
ae t*Rt& m^)u = 0. 


Bessel’s differential equation arises in certain eigenvalue problems and 
{Iml /Ag™x) in = 1,2, ind 


forms a (complete) set of eigenfunctions with weight function x on [0, 1]. Here 2," 
are the positive solutions of 


J,(/2) = 0. 


The Fourier series expansion of an arbitrary function in terms of these eigenfunctions 
is known as a Fourier- Bessel series. 


It was shown how Bessel’s equation arises when separation of variables is applied to 


the membrane (wave) equation and the heat conduction equation in polar co- 
ordinates. 


Incidentally we also included elementary properties of the Gamma function T, and 
we related the Bessel functions of half-integral order to the elementary trigonometric 
functions. 


PDE 14SAQ1-4 


144 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 

Either put m = 0 in Equation (40.5) (M^: page 179): or you can differentiate the series 
expansion for Jo(t) term by term. 

Solution to SAQ 2 


We expand the derivatives on the left-hand side of Equations (40.5) and (40.6) in W, 
so that 


—ntC "7M (o) + (CIO) = tS a(t), 
nt" 7! J,(r) + eit) = td, s. 
Adding t” times the first. equation to t^" times the second, we obtain 
2440) = Jay — Jui, 


whence the result. 


Solution to SAQ 3 


The Gamma function is defined by 


T(m) = [ e "7! dt for m > 0. 
Jo 


- 
« rm-[ etd = [=e = 


o le 


integration by parts gives 


T(2) = f te^! dt = [e] E f Í e^! dt 
Jo o o 


-04r()21; 


and similarly 


r(3) = [ re~ dt = E + af te^! di 
40 o 0 


2042r()22 
(b) T(m + 1) = [ e^" dt 
Jo | 
- len] Í + mf e "n7! dr 
0 o 
= m (m). 


(Note that if m >O then lim e^'t" = 0 as t ~ co, using a result proved in Section 
29.1.3 of Unit M201 29, Laplace Transforms.) 


If m is a positive integer, 
Tim + tl) = mI (n) = m(m — I)F(m — 1) 
= m(m — 1)...2-:1T(1) by induction 
-m(m-1)...2:1-2ml, 
since l'(1) = 1 from (a). 
Solution to SAQ 4 
From the definition of the Gamma function, 


r(4) = [ eiè dt. 


“oO 
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Substituting t = w?, we obtain 


ll 


T) [ e^" u^ Qu du = 2 [ e^" du 
“o Jo 
= 2: irt = m. 
Solution to SAQ 5 
From the series (40.4) in W, we have 


a= v UE SL 
i io 2" * Ekm + k)! 


Differentiating with respect to t, we find that 
d. ce 2k(—1yp*-! 
mit "40-2 — CFC 
al o] P" 2"*?*kYm + k)! 
Ld kam = k 
ss (= Lymn 1+ 2k 
=f Sao 
2, 2m*?57 Ye — T) (im + k)! 


G (— Dnm 120 


ilm m6 2. ntis 
= Sa Mt). 


which is Equation (40.5). 


jim + d +n)! putingk enc 


To obtain (40.6), we multiply the series (40.4) by ¢", obtaining 
n dux 
(= 1m 

Pj = . 

(t) p» 2"* 25k lm + k)! 


We differentiate both sides with respect to t, so that 


dou uno & C2 + gente 
dt Uu) = PN 2"* 3h aam + k)! 
E (= 1k ea 


MELLE Um EET 
ds, - a(t). 


Il 


Solution to SAQ 6 


Let 
u(t) = :f cos(t sin 0) d0. 
T Jo 
Then 
i us f sin 0 sin(t sin 0) dO, 
dt n Jo 
and 
E t a |=- 1 f sin 0 sin(t sin 0) d0 E li sin? 0 cos(t sin 0) d0. 
dii dt X Jo ndo 
Therefore 
E Ir. . . 
i t us (D| + tutt) = -f cos?0 cos(t sin 0) d0 — — f sin 0 sin(t sin 0) d0. 
dii dt 7 Jo T Jo 
Now 


n 


[sn 0 sin(t sin 0) d0 = | -cos 0 sin(t sin nf + Ji cos 0 cos(t sin 0) cos 0 d0 


- f cos?0 cos(t sin 0) dU. 
0 
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Hence u satisfies Bessel’s equation of order zero. Clearly u is bounded as t ~ 0. There- 
fore 


u(t) = AJolt). 


Now 


u(0) = - [ d0 = 1 = Jo(0). 


7 Jo 


so that A = 1 and u(t) = Jo(t), as required. 
Solution to SAQ 7 
We require the\coefficients in the series 
1I-S22YAJQASX  O0<x <1. 
k=l 


We multiply both sides by xJo(,/2,?!x) and integrate over (0, 1). By the orthogonality 
property of the eigenfunctions 


1 
af x[Jo( 2, 0x)? dx = f U = x2)xJo(/2,x) dx, 
0 0 


that is 


MER = [0 - r dx 
0 


using W: page 181, line 17. On the right-hand side the first integral has been worked 
out in the example in Section 14.1.2: 


1 
1 
: Oy) dx = 0) 
Í XJol V2, ®x) dx = i A). 


For the other part, we make the substitution t = /2,!x and obtain 


f JT 
=f X? Jol / A, x) dx = -zol P Jolt) dt 
o 


L^n 


x 


1 A q 
- 2 
a im], t q UO dt * 
1 TE] fi. 
= -70 [eno] -2 f es (t) a} 
^n o Jo 


integrating by parts 


Ant 
i pop mem -i[ eno ar} 
0 


= = rrp AOPE VR) — 2AM}. + 


We have used W: page 180, Equation (40.6) in obtaining the lines marked with an 
asterisk. We now have 


oa VAD 
A/F 


and the solution given in W: page 422 follows. 


Solution to SAQ 8 


We require the coefficients A, in the series 


xi X AJ, GAL IP). 
k=1 
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Multiplying both sides by xJ, (. / 4, ""x) and integrating over (0, 1), we have 
H 1 
A, J, "UA dx = fae Am 
: ))? dx A x" Tlf Ax) dx. 


Using W: page 181, line 17 we find that 


1 
3A, (VAP = | +J (x) dx 


o 
1 NN] 
TIpNWGI +J (0) dt 
[A4 UH +2 f (t) 
1 Jas 
“oe | 


1 asm 


2 
= Te Ine SIP) rue mud 


d 
' e 3: (eet ya (E) at by (40.6) 


] 2 i 
NA Jn LA) = vus Jui A07) by (40.6). 
hi n 
Substitution of A, back into the series gives the solution in W: page 422. 


Solution to SAQ 9 
(a) Equation (40.6) gives us 


q 
BJA = 5 [24]. 


(This expression is valid since J_, is defined by the series for J,, and Equation 
(40.6) was proved for all m in SAQ 5.) Thus 


tJ. a0 = t [ (2) sin | 


Hence 
2\t 
-4U)7 |— im 
J-4(0) B cos 
Alternatively, if you are feeling very energetic, you could show that 


(2k)! 
yk mi 


and hence derive the result by comparing power series. 


kIT(k +4) = 


(b) Equation (40.6) may be expressed as 


my, (= [: 2) [a(o]. 


Applying this formula n times we obtain 


P=, d) = [: L) terso. 


For m = —4 this gives us 
1dy 
= tf- [54 a 
Ji, = 1 |; al it 40] 


Barn 
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Solution to SAQ 10 
We attempt a separable solution 
u(r, 0, t) = RINO T(t). 
Substitution into the differential equation (41.1) gives 
ROT" — e[rors : R'OT + sner| =0. 
After dividing through by ROT we find that 
T" R" R © 
eT Rt iR TRO 
Both sides must be constant, say — À, so that 
T" + CAT =0 (1) 


and 


Rn" ^ R " o " 
R R PO CU 
The second of these equations can again be separated into 
PR" rR e" 
+ ar = - =m’, 
R UR © 


say. (The most convenient form of constant is usually introduced: in this case we 
select m?—a square because we expect to take square roots, with a positive sign 
because we expect the solution for O to be trigonometric rather than exponential.) 
Thus R and O satisfy 


PR" + rR + (Ar? -nm)R-0, (2) 
e" + mO = 0. (3) 


Since the general solution of (1) is A cosc /2t + B sin c / 2t, and 


me, 0,0) = 0. 

we choose T(t) = cos e At. The general solution of (3) is 
Oo(0) = Ao + Bod m=0, 
©,,(0) = A, cos mÜ + B,sinmü — m + 0, 


DL 


and since O(0) must be periodic with period 2m (otherwise it would not be single- 
valued), m can take only the values 0, 1,2,... and By = 0. Finally the solution of (2) 
which vanishes for r = 1 and is bounded at the origin is 


Ju a9). 
the eigenvalues 4 being given by the solutions of the equation 
J4/2) = 0. 


These values are denoted by 2,"". Note that this initial-boundary value problem 
generates a doubly infinite set of eigenvalues. The series (41.3) can now be put together. 


Solution to SAQ 11 
The required solution is given by W: page 182, Equation (41.3). Since 


f cos 0 sin m0 d) = 0 m= 1,2. 


dan 


and 


f cos Ü cos m0 d0 = 0 m= 0,2,3 
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all the coefficients c,,,.d,,, vanish except ¢,,. Thus the solution simplifies to 


Mr 0.0) = Y, c, 0080 J G2, P7) coste / 2,1). 
KSI 


u(r, 0, 1) = c1, cos 0 J,(/ 2, 1r) cos(e 2, 1). 


Finally, by inspection we can see that c, = | in order that the initial condition should 
be satisfied. 


Solution to SAQ 12 


This is the problem of W: page 182 with f(r, 0) = e(1 — r?). The solution may there- 
fore be written in the form (41.3). We require the integrals 


[ dO = 27, 


von 


[ cos m0 d0 = 0 form #0 


“=n 


and 
J sin m0 d0 = 0. 


Thus the solution may be expressed as 


u(r, 0,1) = Y A, Jo / Ax 'r) coste  / 24). 
k=l 


The initial condition gives 


e(l — 7?) = p Ag Jo / 2,7); 


kel 
using the result of SAQ 7 we find that 
IRAN NU EAR 


"0, 1) = 4e a 
u(r, 0,0) eÈ ROUA 


Solution to SAQ 13 


The solution is 


ks — 
u(r. t) = X Au Jo f 2,9 re 7? LR 


n=l 
with initial value 


Es — 


1- = 0) = Y. Ando (A. 
z] 


This is exactly the same problem as SAQ 7 and we can write down the solution as 


JaG/ 2,9) Jo / 2, re hn 
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14.5 APPENDIX 


Equations of Motion for Fluid Flow 


In Section 1.1.3 of Unit 7, The Wave Equation, we derived the equation of motion for 
a compressible fluid moving in one dimension. We now outline briefly the derivation 
of the equations of motion for motion in two or three dimensions. 


First we consider the derivative following the motion. Let the scalar field f represent 
some property of the fluid: we require the rate of change of f(r, t) at a given fluid 
element, whose position changes with time. For a general variation in r = (x, y, z) and 
t we have the first-order Taylor approximation 
ð of àj 
f Y a 


Af S ped aya 
foa eta tx ar 


åt 
à 
= Ar- grad f + Yar: 


thus 


Since we are following the motion we have 


Ar 
lim — = v(r,t 
m Ar (r,t), 
where v represents the fluid velocity. Hence, proceeding to the limit we obtain the 
total derivative following the motion 


which is the three-dimensional analogue of the result obtained in Unit /. 


Next, we consider the generalization of the mass-conservation equation. (This is also 
known as the equation of continuity.) Let V denote the volume contained within a 
simple closed surface S. If AS is a small element of this surface then the volume and 
mass of fluid crossing AS per unit time are respectively vcos 0 AS and pvcos 0 AS, 
where 0 is the angle which the direction of v makes with the normal to AS and p 
denotes the density of the fluid. If we write n to denote the unit vector whose direction 
is that of the normal to AS, then the rate of flow of mass across AS in the sense of n 
is just pv-n AS. Since we adopt the convention in the case of a closed surface S that n 
is always drawn outwards, we find that 


fø *ndS = mass of fluid flowing out of S per unit time 
S 


decrease of mass within S per unit time 


? f: pav 
ôt V 


Now, from the Divergence Theorem (Unit 3, Elliptical and Parabolic Equations), 


ll 


Í pv -ndS =Í div (pv) dV. 
s v 


Consequently, for any volume V, 


i div (py) dV + Í dP aiv =o, 
v y Ot 
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Thus at any point of the fluid we must have 
" dp 
d LEVEL 
iv (pv) 4- ài 0, 


as quoted in Section 14.22, Since 
div (pv) = p div v + v - grad p, 


we can write the expression above in terms of the total derivative as 


dp : 

—+pdivy =0, 

a? 

Finally we consider Newton’s Second Law of Motion for the fluid. Consider a small 
rectangular element of non-viscous fluid with sides of lengths Ax, Ay, Az. The forces 
acting on it are the pressures on its Several faces and possibly a body force F(r, t) per 
unit volume due to gravitation or some other external force. The acceleration pro- 


duced is dv/dt where v(r, t) is the velocity vector r at some point in its interior. The 
x-component of the equation of motion is therefore 


ð di 
F,AxAyAz + pAyAz — [> + x AyAz = T pAxA^yAz, 


+ second order terms, 


where p denotes the pressure. Dividing by AxAyAz and letting Ax, Ay, Az all tend to 
zero, we have 


Combining this with the corresponding equations for y and z we have the vector 
equation 


F d dy 
uem =p— 
grad p = p dr 
In this last equation dv/dt is simply the vector 

dv, dv, dv, 
dt’ dt’ dt}? 


and 
dv, 0v, 
u e" * grad v, + à 
We may condense this as 
dy ðv 
— = (v - grad) yv + —, 
ae cames 


so that the equation of motion becomes 


2 + («grad = — grad p + F. 
t 
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